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1. Introduction. The recent expansion of radio towards ultra short waves has 
aroused a new interest in theoretical problems of electro-magnetism and especially 
in the problem of antenna oscillations and radiation properties. The type of approxi- 
mate discussions used by radio engineers for the case of long wave lengths is of little 
practical value for ultra short waves, where a more rigorous theory is needed, because 
the diameter of the antenna wire can no longer be considered as very small when 
compared to the wave length. 

Some older calculations on rather thick antennas have already been found very 
useful. M. Abraham's! discussion of the vibrations of very long ellipsoids has often 
been referred to. A complete discussion of proper vibrations of ellipsoids of revolution 
may be found in M. Brillouin’s book Propagation de I’ électricitté (Hermann, Paris, 
1904, pp. 314-395) with numerical tables for all eccentricities, from the sphere to 
rather thin ellipsoids. More recently, L. Page and N. I. Adams, and subsequently 
R. M. Ryder? have discussed the free and forced oscillations of all types of prolate 
ellipsoids of revolution; while Barrow,* Schelkunoff,‘ and others have treated the 
problem of the biconical antenna and its free or forced oscillations. Mie and Debye® 
had formerly discussed the free vibrations of the sphere. In most of these papers, the 
theory was based on a computation of the whole field distribution around the antenna 
with the proper boundary conditions on the surface of the antenna. For a perfect 
metal, for instance, the electric field must be orthogonal to the metal surface. 

The aim of the present paper is to emphasize the practical importance of another 
method based on the use of retarded potentials. The principle of the procedure was 
indicated a long time ago,* and the method was recently applied by Hallen and 
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Ronold King’ to the actual computation of antennas. The finite conductivity of a 
real metal can be taken into account, but there are still a few basic questions to be 
discussed, and these will appear more clearly in the problem of a perfect metal with 
infinite conductivity. 

The principle of the method is the following: let us first assume a very thin wire 
and call s a distance measured along the wire. The problem is to find the current 
distribution, Z(s, ¢), along the antenna wire. To such a current, J, there corresponds a 
charge density, o(s, ¢), by the condition of conservation of electricity 


Oa ol 1) 
oF Se Ss \ 
Ot Os 
or, if we assume the following time dependence /(s, t) = I(s)e‘*', 
a Ol 
o(s, d) = — — e*. (2) 
w OSs 


Here, real w means sustained oscillations; while proper oscillations of the antenna 
array will yield complex proper values w, the imaginary part corresponding to radia- 
tive damping. 

An arbitrary current distribution (1), creates an electromagnetic field in the whole 
space which satisfies Maxwell’s equations. This field can be readily computed by the 
method of retarded potentials. In particular, the field on the surface of the metal wire 
can be obtained in this way; and one may then write the necessary boundary condi- 
tion, that this electric field shall be orthogonal to the surface. This yields an integro- 
differential equation which is perfectly rigorous and whose solution is the actual cur- 
rent distribution required. 

Using retarded potentials, one is certain to obtair, at a large distance, a field 
distribution corresponding to a wave spreading out of the antenna. It should be em- 
phasized, however, that the same method can not always be used for the computation 
of oscillations inside a closed tank resonator, where the oscillations are of the type of 
standing waves and have no outside radiation (advanced potentials may sometimes 
be needed too). 

The proper values of this integral equation give the proper frequencies (including 
damping) of the antenna. The same method can be used to study forced vibrations, 
if one assumes an outer electric field acting on the antenna (receiving antenna) or a 
certain electromotive force inserted in the circuit (transmitting antenna). In the sec- 
ond case, one must take into account, for the computation of the retarded potentials, 
the field radiated from a dipole representing the electromotive force. 

Let us discuss the free vibrations of an antenna. The field at a point P is given by 
the well known formulae: 





OV OF , ] o*ds 
A. ss <= « my tee peers V = - . 
Ox ot €o9 o r ; 
s (3) 
oF, oF, ; ¢ itds 
pol - = eee ee SY 46 (p64 - F = pe j : 
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7E. Hallen, Uppsala Univ. Arsskrift 1930, No. 1; Nova Acta, Uppsala, Ser. 4, 11, No. 4 (1938); 
L. V. King, Trans. Roy. Soc. London, 236, 381 (1937); Ronold King and F. G. Blake, Proc. I.R.E., 30, 


335 (1942). 
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h, electric field; H, magnetic field; V, scalar potential; F, vector potential; r distance 
from the element ds on the circuit to the point P where the field is observed; o*, i* 
charge and the current at the time t—r/c; €, uo dielectric constant and permeability 
in vacuum, in non-rationalized units (rationalized units introduce a 1/4 factor in 
the formulae for both potentials). Let us assume an antenna consisting of a straight 
wire along the z axis, extending from z=0 to z=/. We need the z component, h,, of 
the electric field along the wire and must write that this longitudinal component 


vanishes: 








OV OF, w 
-— -_-—-, i* = [(z/)ese-*), kh = — = wren, 
z ot c 
; ba (4) 
- ie** l Ol (z ) e7tkr l e~ikr 
Yu—| — —~s, F, = joe’! f I(2!) —— dz’. 
é9v Jo Oz’ Tr 0 r 


The field at point z is the result of integration over all the points, 2’, of the antenna. 
Finally, we obtain the condition 


= i { OI(z’) d8G(r) 
0 





~ RTC) | dz’ = 0, 


€Qw 02’ zZ (5) 





putting G(r) = , r=|z—2'|. 
This is our fundamental integro-differential equation for the straight antenna. 

One difficulty appears immediately: G is infinite for r=0, z=2’. This means that 
one must take into account the radius of the wire; but when this radius, a, is explicitly 
introduced in the calculation, there is an additional condition to be written for both 
ends of the wire. Here most authors do not attempt to write rigorous conditions; they 
are satisfied with approximations corresponding to the problem of very thin wires. 
They neglect a// but keep terms in Q-', Q-?,- +--+ where 


l 
2 = 2 log —- (6) 
a 


Such a procedure is suggested by the similar approximations used by M. Abraham in 
his discussion of ellipsoids. It should work correctly when 2>14, which means 
l/a > 1000, but could certainly not be relied upon for thicker wires. 

Furthermore, Oseen and Hallen both use the following assumptions: 


0, (2) = 0, current zero at both ends; 


[(z —z')? + a? }*/2, 


A) I(0) 
B) r 


(7) 


The first condition, A, is not quite correct, since there must be a small current at both 
ends in order to charge and discharge the terminal capacities. It is only for the case 
of a hollow cylinder that the current would be exactly zero at both ends; and this 
hollow pipe is a very special case, as shall be seen later. 

The second assumption, B, is explained differently by both writers. Oseen as- 
sumes a current flowing along the axis of the cylindrical wire and computes the field h,, 
Eq. (5), on the surface. Hence his boundary condition (5) is right, but the assump- 
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tion about the axial current is certainly wrong. Indeed, owing to the skin effect, the 
actual electric current, in a perfect conductor, flows along the surface. Oseen assumes 
that the field created by this actual superficial current could be obtained by a ficti- 
tious axial current. This may be right for very thin wires, but the assumption is ob- 
viously wrong for thick wires or for cylinders of large radius. Moreover, it cannot 
be proved that the fictitious axial current satisfies the first assumption A. So Oseen 
hardly justifies the use of both assumptions A and B. 

Hallen takes a different point of view. He starts from the well-known property 
that the current flows along the surface; but instead of computing the field on the 
surface of the same cylinder, he takes the h, field along the axis. This field must 
certainly vanish; and from this fact, Eqs. (5) and (7 B) follow. This necessary condi- 
tion, however, is not sufficient. One may very well have no longitudinal field along 
the axis and still find a longitudinal field on the surface of the cylinder. These approxi- 
mations would probably be all right for very thin wires; but they can certainly not be 
used for thick wires, where B is wrong and A must be replaced by a more elaborate 
condition, in order to take account of the electric currents and charges on the flat 
ends of the cylinder. 

2. Complete statement for a cylindrical wire of finite radius. The antenna is a 
solid cylinder of radius a and height /. The oscillations studied are those of cylindrical 

symmetry where the current ts equally distributed 

around the cylinder and flows along the surface. J(z’, t) 

—1],G) is the total current at 2’, and (1/27) I(z’, t)dg is the 

(i current through a small sector dg (Fig. 1); hence, 

‘«) =L a(z’, t)dz’, Eqs. (1), (2), is the charge per length dz’, all 

around the cylinder, and (1/2m)odz’dgy the charge for a 

small angle dg. For the flat top of the cylinder (z=/), we 

‘ zi.-_\p call J;(p) the total radial current crossing a circle of 

a ae radius p; while o,(p)dp represents the electric charge 

! f between p and p+dp: 

ae 0c; ol, 1 ol; 


; —-—=—, w=——-: (8) 


. ¥ yy , 
oa L f I(z) ot Op w Op 
’ Similar definitions apply for the bottom of the cylinder 











yo. a (s=0) with a current Jo(p) and charge oo(p)dp. The 
( of } ZeO——x positive signs correspond to the directions indicated 
by arrows in Fig. 1. The conditions for continuity of 

| I.) the current around the corners, at z=0 and z=/ read 
af I(a)=-1), — Iola) = 1(0). (9) 
Pia. 1. Let us first study the fields and potentials at a point 


P(z) located on the cylindrical surface. The potentials 
due to currents and charges along the cylinder are the following (e*** factors have 


been dropped): ee ae 
TO ,! ikr 
V.(z) sie +f {= : Ai = de ’ 
r 7 


wf’ S18 Se ge 





(10) 





F.,(z) 
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The currents along the cylinder flow vertically; hence, there are no horizontal com- 
ponents F.,., F., of the vector potential. The distance r is shown in Fig. 1. 

On the flat top of the cylinder, the current flows radially in the horizontal plane; 
hence, the F;, component is zero, but we find horizontal components, F;, and Fyy, 


of the vector potential: 


. a 2r e~tkr dg 
F,,(z) = wo f f Ti(p) cos g — dp, 
0 +o r 2x 
(11) 


. a 2a e~tkr ; dg 
F,,(z) = wo f f Ti(p) sin g — dp = 0, 
0 0 r 2a 
The transverse components F,,, for a point P in the x-z-plane, is obviously zero by 
symmetry: 
, 1 a 2r ol, e7tkr dy 
vic)=— ff 5 Sap (12) 
€pw J 9 0 Op rr 2dr 


and similar formulae for the potentials Fo, and Vo due to currents and charges on the 
bottom of the cylinder. 











The ¢ integrals are of two fundamental types which will now 
be explained in connection with Fig. 2. 








r? = (z — 2’)? + p*? + p’? — 2pp’ cos ¢, (13) 
1 2r e~tkr 
Cinid.t~ he — f dei (14) 
2n 0 T 
1 Q2r e~tkr 
Ob inte f cos ede. (15) 
2a 0 rT 





Fic. 2. 


G, and C; are two functions which will be discussed more fully in section 5. They are 
symmetrical in p, p’ and even functions of z—z’. With these functions, our formulae 


(10, 11, 12) read 


i "wa 
V.(2) ——. — G,(a, Pp; 2 — z’)dz’, 
€qw 0 02’ 


F.,(2) = wo f I(2’)G,(a, p, 2 — 2’)dz’, 
. (16) 


i ° O01; 
Vi(z) = — 


? Gi(p’, Pp, 2— D)dp’, 
€9v Jo Op 





F;,(z) _ wo f I; Ci(0’, — I) dp’, 
0 
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where p=a for the point P(z) on the cylinder. 
We now are in a position to compute the longitudinal field 4,(z) at the point P, 


according to Eqs. (4) and (16). 











h,= - — - — - — - , 
Oz ot 0z 02 
= sls le a,z — 2’) + RI (z’)G,(a, a, 2 — r 
0 02’ Oz . (17) 
*¢ ol; @ 
+f ry = Gi(p’, a, 2 — D)dp’ +f - at Gul a, z)dp’ = 


This is the first integral equation of the problem which corresponds to Eq. (5) for 
the simplified example of a thin wire. It should be noticed immediately that in the 


first integral 


re) 
— G;,(a, 4,2 — 2’) =—=_— — G,(a, 3 — 2’). (17a) 
Zz 02 
This transformation will be very useful, afterwards, in applying integration by parts. 
Another integral equation is obtained by writing the fact that the horizontal field 


component is zero at a point P(p) on the top of the cylinder: 


4 OV. OV, OV 5 OF 1, OF o- re] re 
hp) = -— - — - I vo T seen 4 
Ox Ox Ox ot ot Ox dp 
€ow : a 2 
er ee h.(p) = f ae —G;,(a, we l—2z "\d2’ 
2 0 02’ Op ; 
(18) 
ol; @ 
+f — 5p ales os 0)dp’ +f°= — * Gils, 1) dp’ 
0 Op’ dp 


Pe f [T1(0")Cu(p", p, 0) + To(o’)Cu(0’, p, 1) Jdp’ = 0 
0 


A similar equation could be written for the bottom of the cylinder; but this is actually 


not needed, since it reduces to (18) by reason of symmetry. 
The proper oscillations of the cylinder can be divided into two groups: 


symmetrical oscillations I;(p’) = Io(p’), I(l—z) = —I(z), 
ol(i—z) 9al(z) 





’ 





Oz Oz 
(19) 
antisymmetrical oscillations I:(p’) = —Io(p’), J(/—z) =1(z), 
ol(l — z) dl (z) 
dz 7 0z 


These two types will be discussed together in the following formulae. The upper sign 
corresponds to symmetrical and the lower sign to antisymmetrical vibrations. 


1943] THE ANTENNA PROBLEM 207 


3. Discussion of the first integral equation (17). Wave propagation along the 


cylinder. Equation (17) can now be written in the following way: 
. -G,(a, a, 2 — 2’) + RI (2’)G,(a, a, 2 — “| dz’ 
“ie ol, oa : ’ ? 
i —_— IGi(p, a, 2 —1) +G,(p’, a, z) |dp = R(Ji, z). (20) 
Op OZ 


[he left hand integral contains only vertical currents, J(z’), along the cylindrical 
boundary; while the right hand terms, R, show the coupling between these vertical 
ind the currents or charges on both flat ends of the cylinder. 


rate the left hand integral by parts, starting from 0G;/02’: 





eirg2] 7 aT z’=l 
| |\— + | 1|G.(a, a,z — 2')dz’ = ——G,(a, a, — 2’) + R(Ii, 2) 
Jo LOZ * a Oz z'=0 
‘... 
(—) G.(c,a,2—D ¥G,(a, a, 2)] + RU, 2). (21) 
\ oz JI tus 
[his new formula has been obtained without any approximations. Let us now make a 
simpl ptions, in order to get a better understanding of the meaning of 
Ss ed TI" 
For a very thin and long wire, 1 >a, we may neglect the R(J;,2) term, as both charges 
and currents on the flat terminals become very small. Furthermore, at a certain dis- 
nce fron terminals, G,(a, a, z—1) and G;,(a, a, 2) are also very small, since G; 


ximately like 1/r for large distances. The only important term is the 


gecreases app! 


one on the left, which has the obvious solution 
a7] : w 2a 
— + kJ (2’) = 0, k=-—=-—- (22) 
Oz * Cc r 


This shows wav 
wire. This result is obtained under the assumption [><a and without any restriction 
about the wave length \, which can ‘be of the order a or even smaller; but it holds 
liurn part of the wire, far away from both ends.® 

This shows the connection with the usual elementary theory of antennas. The 
classical discussion® starts from the assumption of sinusoidal standing waves along 
the wire, which cancels out completely the left hand integral in equation (21). Then, 
using this current distribution, the longitudinal field along the wire may be computed; 


and according to (17) and (21) it comes out as 


h,(z) = J {(=) IG.(a, a, 2 — 1) $ G,(a, a, z) | + R(T, at. (23) 


ve propagation with the velocity of light along the major part of the 


only for the me 


€n® 02’ 2’=! 

8 It should be emphasized, here, that our discussion is limited to the case of oscillations with cylindri- 
cal symmetry (see beginning of Section 2). Vibrations with nodal lines parallel to the axis are not included. 

*L. Brillouin, Radio-électricité, loc. cit. 

J. A. Stratton, Electromagnetic theory, McGraw-Hill, New York, 1941, pp. 455-460. Stratton uses 
rational units, hence a 1/42 factor before the integrals, and he uses the opposite sign before 4. 
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This plays the role of a small additional average impedance Z along the antenna, 
which can be defined by 


Zr? -{ h,(z)I(2)dz. (23a) 
0 


The real part of Z is called the radiation resistance, Z,, and the expression Z,]? repre- 
sents the energy, W, radiated at large distance (see Stratton, p. 458), from which the 
damping of the antenna oscillations may be computed. For a very thin wire, one may 
neglect the term R(J;, z), which represents the role played by the currents and charges 
on both flat terminals of the wire; and one may take for G, the expression (1/r)e~**" 
as in Eq. (5). With these approximations, our equation (23) becomes identical with 
Stratton’s Eq. (76a), p. 457. 

It should be noticed that Eq. (23) is physically wrong, as we know in advance 
that the longitudinal electric field along the wire is zero. These equations: (23) and 
(23a) merely represent a second approximation in a system of successive approxima- 
tions starting from (22). An attempt will be made, in the next section, to build up a 
consistent system of approximations of similar structure. 

Returning now to Eq. (20), we may try another integration by parts, starting from 
0I/dz’, which yields 


: : 0°G; : 
f T(z’) |< + kG,(a, a,z — “| dz’ 
0 OZ ° 
}z/a=l 


re) 
= T(z’) 5,7 Galas a,i= 2’) + R(Ii, z). (24) 
oz 


Z lz’=0 


Let us again discuss this equation for a very thin wire. The term R(J;, 2) represents 
the role of both terminals and may be neglected, J(z’) is zero at both ends (z’=0, 
z’=I1), and consequently all the right hand terms are zero. This transformation is 
very closely connected with the one used by Schelkunoff and Feldman” in a recent 
paper. These authors discuss the problem of forced vibrations in a transmission an- 
tenna, instead of the free vibrations which we have in mind. They use both approxi- 
mations (7A) and (7B) of Oseen and Hallen and take for G the simplified expression 
(1/r)e-**", Eq. (5). These approximations may apply for a very thin wire. Further- 
more, they split the (1/r)e-** function into its real and imaginary parts before per- 
forming the integration by parts. Their final result is actually identical with the one 
derived from the elementary theory and Eq. (23). This is not surprising, as both 
methods are very closely connected. 

4. Principle of a method of successive approximations. As stated in the preceding 
section, it seems possible to build up a method of successive approximations in order 
to solve Eq. (21) along a way rather similar to the one followed in the classical ele- 
mentary discussion. 

First of all, we may split the integro-differential equation (21) into an integral 
equation and a differential equation, by writing: 


10S, A. Schelkunoff and C. B. Feldman, Proc. I.R.E., 30, 511 (1942). 
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l 
f F(z')G,(a, a, 2 — 2’)dz’ = R'(z), (25) 
0 
where R’(z) = (01 /02’).-:[Gi(a, a, z—l) $G;(a, a, 2)]+R(Ti, 2), 
o7y aie oe 
32/2 + RI (z ) = F(z ). (26) 
02°” 


The first equation is an integral equation of the first kind, with the kernel G;(z—2’). 
Its solution can be written with the help of the resolving kernel H;(2’—2”), which 
satisfies the following conditions 


l 
f G.(z — 2')H,(2’ — 2”)dz’ = 5(z — 2”), (27) 
0 


l 
F(z’) -{ R'(2")H,(2’ — 2")dz", (28) 
0 


where 6 means a delta function. Hence, the first question is to build up the resolving 
kernel H;, a problem for which some general methods have been developed. This 
being done, we are left with Eq. (26) to which we apply the usual Rayleigh-Schréd- 
inger method of successive approximations. Let us first notice that the G, function 
becomes very large for z=z’ which, according to (27), means that H; is small. Thus 
we may rewrite (26) and state explicitly by an ¢€ coefficient the smallness of the right 
hand term: 


49 


07] 
aa t HIG) = 92), F =o. (26a) 
Oz“ 


Then we use the following expansions: 








T(z!) = Io(z") + €1y(2") + €2(2’) +++, 
. 2 2 (29) 
k? = ko + exi te x2°°? 
and obtain the successive approximations: 
077 2 
— + kolo = 0, 
02°? 
oF, 
a RTH (30) 
ae 
[Arms + Role =_— xolo _ xii “2 
Oz"? 


Io is a sinusoidal function, as in the elementary treatment, 
Io = A sin ko(z’ + £) 


where the ¢ constant is necessary in order to give a small but finite value for the 
current J» at the bottom of the cylinder (z’ =0). This is needed for the junction with 
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the currents on the lower flat end of the cylinder. By symmetry, the correction at the 
upper end must also be ¢; hence, 
é Ao 2 
Ro(l +- 2¢) = nt, 1+ 2¢ = >? Ro =—-e (31) 
é 0 
The constant ¢ will be determined by means of the second integral equation (18) 
for the flat terminals. Now let us turn to the second equation (30). As is well known, 
it is necessary for the right hand term to be orthogonal to the solution of the homo- 


geneous equation, which means 


I 
f sin ko(z’ + §)[- xilo + ¢|dz’ = 0 
0 


or 
l l 
Axi f sin? ko(z’ + f)d2’ = f yg sin ko(2’ + £)dz’. (32) 
0 0 


This yields the correction x: to the proper value &. It is readily seen that equation 
(32) is very similar to the relation (23a) used in the elementary theory to obtain the 
average “radiation resistance” of the antenna and thence the damping coefficient in 
the proper oscillations. The important point, however, is that equation (32) contains 
y, which is not R’ but is computed from R’ by means of (28)—(26a). 

Once x: is obtained, the second equation (30) can be solved; then x2 is first com- 
puted by a similar orthogonality condition, and so on. Hence, the whole procedure 
should yield a solution along lines parallel to the elementary treatment and show 
how far the usual formulae can be trusted. 

We may already go one step further and write the general expression of the func- 
tion F(z’) on the basis of Eqs. (25) and (27): 


( 


(a!) 


Il 


I 
ey F(z’) = f R’(2")H;(2’ — 2”)dz" 
0 
Oly a jee 
- (=) [6(2’ — 1) ¥ 6(z’)| +f R(1i, 2”) H,(2’ — 2”)dz". (33) 
OZ z’=l 0 


The 6 functions appear here automatically, because G; is an even function of (z—2’), 
and so is Hi; for z’—2”; hence, the integral in (28) comes out as 


feu — I) His — 2")dz" = f exc — 2")H,(2" — 2')dz" = 8(z’ — 1) 


x 


according to (27). 
We can use the new expression (28) for the discussion of some simplified examples. 


Let us start with the wire of vanishing radius. The whole R(J;, 2”) term, which repre- 
sents the terminal effect, drops out; and we are left with an equation 


a7] cas al 
— + (2!) = Fe’) = (—) [i(z’ — 1) F &(z’)). (34) 
z’=l 


4./ 
Oz 


o 
2 


from (26) and (28). The condition on both terminals is obviously J(0)=J(/) =0; 
hence ¢ =0 in (31), which results in the following equation: 
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Io = A sin hoz’, ky = nr/I. (35) 


n =2m-+1: symmetrical oscillation, sign — =(—1)” in bracket, 
n =2m: antisymmetrical oscillation, sign-++- =(—1)”" in bracket. 


l 
xi — = bv f [(— 1)"8(2’ — 1) + 6(z’)| sin kos’dz’ 
0 
= ko|(— 1)" sin (Rol) + sin (0) |] = 0 (36) 


which gives no damping at all. The physical explanation is the following: a finite 
amount of energy is radiated per second; but this does not mean any damping of the 
oscillations, because the energy accumulated in the field around the wire is infinite. 
As a matter of fact, both electric and magnetic fields are infinite as 1/r near the wire 
of infinitely small radius. The square of the field is of the order 1/r?; and the energy is 
obtained by multiplying by 2zr dr and integrating with respect to r, which gives 
logarithmic infinite terms. The situation is similar to the one obtained in a circuit 
with infinite L, zero capacity, and finite resistance R, which yields a negligible damp- 
ing coefficient R/2L. 

This shows the difficulties involved in the assumption (7A), as put forth by Oseen 
and Hallen. When such a condition is used in the rigorous Eqs. (25), (26), it leads 
directly to (36) and yields practically no damping. 

Such is also the case for a hollow cylinder. Here again, there is no end effect, no 
terminals, no R term, and condition (7A) holds good. The whole procedure from (34) 
to (36) repeats itself and shows again no damping. Of course, the G, and H; functions 
would differ materially in both cases; but these 


functions have been eliminated from Eq. (34) and ais ELECTRIC 
ye SY FIELD 


finally drop out. 
The explanation is similar to the one given for os i 
the thin wire, but not quite so obvious. The prob- I 1: I: i 
lem of a hollow cylinder of indefinitely small thick- e| | ° | - 
ness must be considered as the limit of a cylinder CURRENTS 
of finite wall thickness, as represented in Fig. 3. 
On such a cylinder, one should take into account, 
separately, a current J; flowing along the external 














surface of the cylinder and another current J; along ‘ i 
the internal surface. At the limit, these two cur- E , 
i a 
rents merge into a single one, for which the theory 
Fic. 3 


indicates a sinusoidal distribution. Hence, for a 
cylinder of finite thickness, there certainly is a 
current flowing around the edge of the cylinder, as shown in Fig. 3. On this edge, 
one must also consider the electric charge; and this results in an accumulation of 
electric fields and of electro-magnetic energy near the cylinder, while the energy radi- 
ated per second at large distance remains finite. Hence the damping becomes negligible. 

The result is general and applies for any hollow cylinder of indefinitely small thick- 
ness, whatever the shape of the cross-section might be. The field distribution inside 
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the cylinder should correspond to a superposition of Ey waves (transverse magnetic) 
and should show a strong decay from both ends down to the middle part of the cylinder, 
especially when the diameter of the cylinder is small compared to the wave length. 
These two simple examples show the importance of the role played by the shape of 
both terminals and the danger of using assumptions like (7A) or (7B). 
5. Some important formulae. We have introduced in (14), (15) two fundamental 








functions: 

“ee P 1 ae g thr , 

Gx(p, p, $) = — dg, f=s-—s', 
2r 0 r 
1 2r e~ikr 

Ci(p, p’, ) = —{ cos gdg, (37) 
2a 0 Tr 

r? = ¢2 +4 p? + p!? — 2p’ cos y = g — 2p cos g, q = 2p, 


Ul 


q =? + p? + pp”, Pp = pp. 


From these relations, we see that G,; and C; depend upon p, p’, ¢ only through the 
two combinations p and g. Furthermore, it is easily proved that 


OG; . OC, 1 27 1 @ (<< 





— 


) cos gdg. (38) 


Op 0g 2rJo «6Urd«€OOrN\ hr 


C; and G; being both zero at infinity, this can be written as 


1 " OG, 
Cs = —{ weet dq. (39) 
2J, op 


These integrals are closely connected with the complete elliptic integrals K and D,” 
as is seen for a thin wire when the radius a is small compared with the wave length 
(ka small). The following expansions can be used: 


r=V/q— 2pcosy = Vg + [V9 — 2p cosy — Vq] 

| agi Id Bal a (40) 
e~ikr = e-ikvat1 — ik[y/qg — 2p cosy — Vq]--:: }. 

The bracket [ | is of the order of magnitude of a, and its product when multiplied 

by k is small: 


. 17% 1+ ik/q 
G, = eikva — f _ wade — ik} 
2nJo [gq — 2p cos g]?/2 


17% 1+ ikv/g 
C, = e«** i— f — : cos ygdg —0+-:---: \. (41) 
2eJo [gq — 2p cos ¢]}/? 








We may write 
g— 2pcos¢ = (q + 2p)(1 — x? sin? y) 
4p g-T 
=-=——, y=. (42) 
q+ 2p 2 


K2 


Hence 


1! E, Jahnke and F. Emde, Tables of functions, 2nd ed., Springer, Berlin, 1933, pp. 127-145. 
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i deg 2 er dy _ 4 K(.) 
> [g—2pcose]® [gt 2p]}?I_s2 [1—asinty]? [g+ 2p}? 
and 
[2 1+ ikV/q 
G, = etka] —. ‘+ ove ikv/q K(x) — ik--> |. (43) 
m [q+ 2p]? 


When £->0, the variables g and p retain finite values; but when at the same time 
p=p', g=2p, then x is 1 and K is logarithmically infinite. This could easily be fore- 
seen and does not make any special trouble in the integrations. The second integral C; 
is transformed in a similar way: 


- cos gdy 2 {~ 2 sin? y S| a 
0 [g—2pcosy]!? [g+ 2p]? J _./2 [1 — x? sin? py]? 


[2D(x) — K(«)], (44) 








4 

[g + 2p]? 
. {2 1+ iky/q 

Cr = e7tk¥a ¢ —.—_____—_ 
w [q+ 2p]! 


These approximate formulae should be used for a thin wire and represent the first two 
terms in an expansion when a/A is small but not negligible. For the fundamental 
vibration, \ is of the order of 2/ (twice the length of the antenna). Hence using the 
expansions (43), (44), one should be able to go one step further than Oseen or Hallen, 
who completely neglected a// and were satisfied with keeping terms in Q-', Q-?, where 


[2D(x) — K(x)]--- \. 


l 
Q = 2 log —- (6a) 
a 


This parameter comes in, when integrations are performed on D and K for « near 1, 


small; ¢=s-— 32’. 


_—g-% +(o-P')? 

qt+2p + (p+ 9’)? 
This happens when z and 2’ are nearly equal for two points on the cylindrical surface 
p =p’ =a. It happens again for two points on one of the flat terminals, when z=2’ =0 
or ], and p is nearly p’. In such cases, K and D are represented by the following ex- 
pansions (Jahnke-Emde, p. 145) 


x2? = 1 — x? 








’ 


K=A+——«?+--, D=A-142(A-Dx'?---, 





me om (45) 
(s — 2’)? + (p — 9’) 
@— 2)? + (+0) 
Integration and averaging process carried out on A will introduce the parameter Q. 
Finally, let us discuss the dependence on k of both functions G, and C;. From 


the definition itself (37), it is seen that both functions can be expressed in terms of 
G;, C, corresponding to k=1, 





4 
A = log— = log 4 — 3 log x’? = log 4 — 3 log 
kK 
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k e ikr 

Gi = -f — dg = kG,(kr) hence: 
2a kr 

Gi(p, p’, $) = kGilkp, kp’, kf) = kGi(k?q, kp), 

Cr(p, p’, $) = RCi(kp, kp’, RS) = kCi(R*q, B*9). 


(46) 


The same decomposition can be seen from the expansions (41). 

6. Conclusions. The preceding sections show clearly the importance of the role 
played by both end-surfaces, whose exact shape should be taken into consideration 
very carefully. We have shown, on the example of plane terminations, that the prob- 
lem consists in finding ‘two unknown current distributions, one for the cylindrical 
surface and one for the (symmetrical) terminal surfaces, and this requires solving 
two integral equations. This is the essential difference from the problems of the proper 
oscillations of one closed algebraic surface, such as an ellipsoid. For plane termina- 
tions, a complete study of equations (17) and (18) should be affected, and the suc- 
cessive approximations should be worked out simultaneously on both equations. 
Other shapes of end-surfaces, like half spherical or half ellipsoidal terminals, would 
certainly yield quite different results. A discussion of this problem is not attempted in 
the present paper, the aim of which was merely to offer a precise statement of the 
mathematical theory of antennas and to emphasize some difficulties which seemed 
to have been overlooked by previous authors. 
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STABILITY OF COLUMNS AND STRINGS UNDER PERIODICALLY 
VARYING FORCES* 


BY 


S. LUBKIN anp J. J. STOKER 
New York University 


1. Introduction. It is a well known fact that a rigid body hinged at one end and 
standing vertically can be put into stable equilibrium by applying a vertical periodic 
force of proper frequency and amplitude at the lower end. The differential equation 
for small oscillations of the rod is a linear homogeneous equation with a periodic co- 
efficient—it is a Mathieu equation if the applied force is a simple sine or cosine func- 
tion of the time. Stability of the rod would require that all solutions of this equation 
be bounded; it is found that this is the case if the frequency and amplitude of the 
applied force are properly chosen. A more complicated problem of the same general 
type in a system with more than one degree of freedom has been considered by 
G. Hamel [4]'; linear differential equations with periodic coefficients play the essen- 
tial role in this case also. 

We shall be interested here in analogous problems in elastic systems with infinitely 
many degrees of freedom. One of these is the problem of the column under periodic 
compressive force F(t) applied at the ends of the column.? The analogue of the 
problems mentioned above would be as follows: the force F(t) consists of a constant 
part P plus a periodic part H cos wt. Suppose that P were a compressive force larger 
than the lowest compressive load (the Euler load) for which the column in the 
original unbent position is instable. The question is, then, whether or not H and w 
can be chosen in such a way that small motions in the neighborhood of the unde- 
flected position are stable ones. We shall see that this can always be done, though, 
as one would expect, the quantity H must be chosen so that the total force F(t) falls 
below the Euler value during at least part of the time. However, the time average 
of F (over a cycle) may be very much larger than the Euler load. On the other hand, 
it is quite possible that the column may be instable when P is a compressive force 
smaller than the Euler load or when P is a tension rather than a compression, if 1 
and w are properly chosen.’ From the point of view of the practical applications these 
latter possibilities are certainly the more important ones. For the case of the column 
with pinned ends we give diagrams which make it possible to decide whether the 
column is stable or not under any of these circumstances. The stability of the 
stretched string under a tension which varies periodically in time is also considered. 

In all of these problems the Mathieu equation’ (more properly, a sequence of 


* Received April 9, 1943. 

1 Numbers in square brackets refer to the bibliography at the end. 

2 A special case of this problem has been treated by I. Utida and K. Sezawa [16]. 

’ Analogous problems for plates under loads in the plane of the plate have been considered by 
R. Einaudi [1]. 

* We consider always that the applied forces are simple harmonic functions of the time—otherwise 
we should have to deal with the more general Hill’s equation. 
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Mathieu equations in the continuous systems) plays a central réle, since the decision 
as to stability depends upon the character of the solutions of such equations. For this 
reason a brief summary of the main facts concerning the solutions of the Mathieu 
equation is included here. A brief treatment of the Mathieu equation with a viscous 
damping term added is also included because of its importance for the stability 
problem. 

2. The column under periodic axial forces at its ends. We make the assumptions 
that are customary in dealing with the transverse oscillations of thin rods. Of these, 
the principal ones are: 1) the rod is an initially straight uniform cylinder, 2) the 
lateral deflection w (Fig. 1) and the cross sectional dimensions of the beam are small 
in comparison with the length /, 3) all stresses remain below the proportional limit, 


F(t) F(t) 
a {wane eS = i See See ee Se ee ee eS ee *, }———+ 











4) the effects of shear and rotary inertia are negligible.5 In addition, we assume that 
the column is subjected to axial forces F depending on the time ¢ and applied at 
the ends of the column; these forces are counted positive when they are tensions. 
With these assumptions the differential equation for the lateral deflection w(x, f) is 
well known to be as follows: 

04w 0*w 0*-w 


EI - eae Fi) —— + w-— = G, (2.4) 


Ox’ Ox? ot” 


In this equation E and J are Young’s modulus of the column and the moment of 
inertia of its cross section, and m is the mass per unit length. In what follows we 
assume always that F(t) is given by 


F(t) = P + H cos 2rft; (2.2) 


i.e., it consists of a constant part plus a harmonic component of amplitude H and 
frequency f. 

It should be pointed out that the derivation of (2.1) involved a tacit assumption 
not included among those enumerated above. This was that the forces F(t) applied 
at the ends of the column result in forces throughout the column which are, to a 
sufficiently close approximation, independent of x. We proceed to show that this 
assumption is warranted under the circumstances normally encountered in practice. 
The differential equation for the longitudinal displacement u(x, t) of the rod is 

0°u O*u 
re Sp 
Ox* ot 


(2.3) 


in which p is the density of the rod. The total force F transmitted through any cross 


section of the rod of area A is given by 


5 These effects could be taken into account without difficulty, but nothing new in principle would 


result. 
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Ou 
F =AE—.- (2.4) 
Ox 
We assume as boundary conditions 
u=0 at x«=0, (2.5) 
and 
Ou 
F=AE - = P+Hcos2xft at x =/1/2, (2.6) 
x 


the origin of coordinates being taken at the midpoint of the red in order to take 
advantage of symmetry. We seek the forced oscillation and neglect the free oscilla- 
tion. The result for the quantity F is readily found to be 


F(x,t) = P+H tame 2nft (2.7) 
— cos (Al/2) adie 
with 

A = 2nf(p/E)'*. (2.8) 


It is convenient to introduce the fundamental frequency fo of the free longitudinal 
vibration of the rod which has a single node at the center. This is given by 


fo = (1/21)(E/p)*"”. (2.9) 


Upon introducing this into (2.7) we obtain 

cos (rfx/fol) 
————- co 
cos (rf/2fo) 


If f is small compared with fo it is clear that F will be nearly independent of x. For 
steel or aluminum (E/p)'/?=17000 ft./sec., while for brass, concrete, stone, or wood 
this quantity is about 12000 ft./sec. For any column of usual length fo will therefore 
be of the order of 500 cycles/sec. or more. Hence if the applied axial force F(t) is one 
of frequency below say 50 cycles/sec. it is reasonable to assume that the variation of 
the axial force with x may be neglected. 

We introduce new independent variables replacing ¢ and x in (2.1) by the equations 


F(x, =P+H s 2nft. (2.10) 


&= 2rft and §£ = xrx/l. (2.11) 

In addition, it is convenient to introduce new parameters as follows: 
Pg = vw’ El /I’, é9 = Ps/EA, (2.12) 
p = P/Pz, h = H/Pz. (2.13) 


The quantity Pg is the negative of the Euler load for the column and € is the tensile 
strain due to that load. The quantities p and h are the ratios of the constant part and 
of the amplitude of the oscillating part of the applied load to the negative Euler 
load. With these new quantities the differential equation (2.1) becomes 


O4w 0*w O*w 
— — (p+ h cos #) e + (f7/fi 0) 


oe 
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The quantity fp is the fundamental frequency of longitudinal vibration of the column 
given by (2.9). 

The general problem which we wish to investigate can now be stated: for given 
boundary conditions there are certain values of p, h, and f for which all solutions 
w(£,#) of (2.14) remain bounded when arbitrary initial conditions are prescribed and 
other values of these quantities for which unbounded solutions exist. In the former 
case we say that the column is stable and refer to p, h, and f in this case as stable 
values. Our problem is to separate the stable from the instable values of p, h, and f. 

We do not solve the problem in this generality; we choose rather a special case 
with regard to the boundary conditions to be imposed. 

3. Formulation of the stability problem for the column with pinned ends. The 
boundary conditions we choose are those corresponding to the case of a column with 
pinned ends; that is, we assume that the deflection w and bending moment 
M =EI(0’w/dx*) are both zero at x=0 and x=/. We have, therefore, as boundary 
conditions for (2.14): 

0?w 
w=—=0 for §=0 and =z. (3.1) 
0g? 

These boundary conditions can be satisfied by taking for w a solution in the 


form of a Fourier sine series: 


w= > F,(8) sin né. (3.2) 
n=1 
The series (assuming that it converges properly) is a solution of (2.14) provided that 
the function F,() satisfies the differential equation 


+ (au + Bu C08 PFs = 0, rong ee eae (3.3) 
in which 
an = n° fi/€o/f?)(n® + p) (3.4) 
and 
B, = n*(f? €o/f*)(h). (3.5) 


The quantities f, fo, €0, , and # have been defined by equations (2.2), (2.9), (2.12), 
and (2.13) respectively. The differential equation (3.3) is, of course, a Mathieu 
equation. 

We can now see why the choice of the boundary conditions (3.1) brings with it 
essential simplifications. To begin with, it is not possible to separate the variables in 
(2.14) in the usual way: if we insert for w in (2.14) an expression of the form 
w=f(£)F@) we do not obtain a pair of ordinary differential equations for f and F 
alone. By assuming for w the special form given in (3.2) we are able to satisfy (2.14) 
by virtue of the fact that only even ordered derivatives of w with respect to £ occur 
in it. This form of solution is, however, not useful for boundary conditions other than 
those given by (3.1). The reason for this is as follows: since w satisfies (2.14) we 





6 The problem can be solved for other boundary conditions, but with much more difficulty. It is not 
possible, for example, to make use of the theory of the Mathieu equation in other cases. For a possible 
approach, see R. Einaudi [1], and S. Lubkin [8]. 
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must require that 0‘w/d& be continuous, since w and 0°w/d& (the bending moment 
within a constant factor) should be assumed continuous on physical grounds. But the 
sine series (3.2) can be differentiated four times with respect to & if, and only if, 
wand 0’w/d& vanish at £=0 and =7, the end points of the column.’ 

Our definition of stability requires that w(t, 3) be bounded for 0S0< @ when 
arbitrary initial conditions are prescribed. Hence we must require for stability that 
all solutions F,(@) of (3.3) for n=1, 2,3, --+- and OS0< © remain bounded when 
arbitrary initial conditions are prescribed. This is, of course, only a necessary condi- 
tion for stability. However, we show in an appendix that the Fourier series (3.2) 
will, roughly speaking, converge for all 3 if it converges for } =0 and if each F,@) 
is a stable solution of the Mathieu equation. Such a question does not arise in the 
more usual type of initial value problem, since the functions analogous to F,@) 
are generally of the form e™°(A,, cos n0 +B, sin nv), r, =0. 

4. The Mathieu equation. The problem of the stability of the column with pinned 
ends has been reduced to that of determining whether all solutions of the Mathieu 
equation 


wi + (a + B cos v)F = Q, (4. 1) 
dv? 

i.e., of Eq. (3.3) without subscripts, are bounded for given values of @ and 8 or not. 
We summarize briefly the known theory of this equation in so far as it is needed 

for our purposes; more extended discussions and proofs can be found in the pamphlets 
of M. J. O. Strutt [15] and P. Humbert [5], and in the books of E. L. Ince [7] and 
Whittaker and Watson [17]. We have also made use of papers of S. Goldstein [2], 
E. L. Ince [6], and M. J. O. Strutt [14]. The notation we have chosen for the 
Mathieu equation has been taken to fit our problem; we compare it with the notation 


used by others: 


| Ince and 








S t xoldstei wae , re 
— Goldstein Whittaker and Watson Here 
u y y F 
2x 2x 2x ov 
/4 a a/4 a 

—h?/2 —4q 4q 8 


It can be shown (theorem of Floquet) that there exist in general two linearly 
independent solutions F; and F, of (4.1) which satisfy the relations 


Fi “fe 2) = K,F (8), 
Fd + 27) = K.F (8). 


(4.2) 


The quantities K; and Ke are either conjugate complex or real constants which 


satisfy the relation 


7 The analogous problem of the rectangular plate with simply supported edges can be treated in the 
same way as the column with pinned ends. The only essential difference would be that the relations corre- 
sponding to (3.4) and (3.5) would contain more parameters. 
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Hence all solutions of (4.1) will be bounded only if 
| Ki| =| K:| = 1. (4.4) 


In case (4.4) is not satisfied, it follows from (4.3) that K; and Kz are both real—a 
fact of which we make use later on. For certain values of a and Bf there exist solutions 
for which the values of K are +1 or —1; such solutions are therefore periodic of 
period 27 or 4m respectively.’ The pairs of values (a, 8) for which such periodic 
solutions of (4.1) exist can be shown to fill out curves in an a, B-plane which divide 
that plane into “stable” regions in which (4.4) holds and “instable” regions in which 
it does not hold. The boundary curves themselves belong to the instable region, the 
general solution of (4.1) corresponding to (a, 8) on such a curve consisting of the 
sum of a periodic function plus # times a periodic function. Fig. 2 indicates these 
regions, the stable ones being shaded. 
It is of some interest to note that in the stable regions the relation 


a+ |B|>0 (4.5) 


must hold since otherwise d?F/dd*? would always have the sign of F and a solution 
not identically zero could not remain bounded for 3+ © as well as for J-— ~; 
this would mean instability since F(—d) is evidently a solution of (4.1) if F@) is. 

The stable regions are connected at the points a=k?/4, B=0, k=1, 2, 3,---, 
for which the solutions of (4.1) are evidently bounded. As indicated earlier, the 
boundary curves separating stable and instable regions are characterized by the fact 
that a periodic solution of period 27 or 47 exists for any pair of values (a, 8) on such 
a curve. This can be made the basis of a method (due to Ince [6]) for determining 
these curves, as follows: a Fourier series with undetermined coefficients is assumed 
as a solution of (4.1). Upon substitution in (4.1) an infinite set of linear equations in 
the coefficients is obtained, each of which involves only three successive coefficients. 
Each equation may then be solved for the ratio of two successive coefficients in terms 
of the next higher or of the next lower coefficients. By successive substitution in 
these relations one is in this way led to two expressions for any such ratio, one of 
which is a finite and the other an infinite continued fraction. By equating the two, 
a relation between a and @ is obtained which holds at the boundary points separating 
the stable and instable regions. For a given value of 6 and with @ ranging from 
—« to + one comes first upon the boundary curve Cy which begins at a=0, 
B=0 (cf. Fig. 2)®; the periodic solutions corresponding to points on this curve are of 
period 27. Following this, the next two curves, C; and Si, starting at a=1/4, B=0 
correspond to solutions of period 47, followed by two, S2 and C2, starting at a=1, 
8 =0 corresponding to solutions of period 27, etc. The letters Cand S refer to develop- 
ments in cosine series (for the even solutions) and in sine series (for the odd solutions). 
The points between two successive curves for which the periods of the corresponding 
solutions are different are stable points. For small 8 the boundary curves are given 
by the following expressions, solutions of type C2, and S2, having the period 27, 
while those of type C241, S2e41 have period 47: 

® For a given value of 8, say, the problem of determining values of a for which such solutions exist is 


obviously a linear eigenvalue problem. 
* Essentially the same figure appears in the book of Strutt [15]. 
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Curves of type C2; and S2, have contact of order 2k at the points a=k*, 8=0, while 
curves of type C2x41 and S2x41 have contact of order 2k+1 at the points a = (2k+1)?/4, 
8=0. This behavior is clearly indicated in Fig. 2. A table of values of a and 8 for 
points on these curves is given at the end of the paper. These values were calculated 
by means of the procedure outlined above and were checked against values given by 
S. Goldstein [2] and E. L. Ince [6] where possible. 

For large positive values of a the points for which |8| <a are stable except for 
very narrow strips which lie near the lines a=k?/4. For large values of 8 it has been 
shown that all boundary curves tend to have the slope —1 (for 6>0). Thestable 
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regions are in general very narrow for a<0 and grow narrower as | | increases. 
These observations are all borne out by Fig. 2. 

5. The stability of the column with pinned ends. We may now conclude that the 
column with pinned ends will be stable only if the applied force F=P+H cos wt 
is such that all points (a@,, 8,) given by (3.4) and (3.5) fall within the shaded region 
of Fig. 2. In other words, a set of values (p, h, f) is stable only if every point of the 
sequence (a,, 8,) determined by (p, h, f) is stable. 

Suppose, for example, that P=Pz, (i.e., the steady part of the load is a tension 
equal in value to that of the Euler load) and that the harmonic part of the load has 
a frequency f=fo(€9/2)'/*. We find that a,=1 and that the column (it is, rather, a 
tensile member in this case) is instable even for small amplitudes H of the oscillatory 
part of the load (i.e., for | B1| small), since the points (1, 6:1) are clearly seen with 
reference to Fig. 2 to be instable if | 6; | 0 is small. We could expect the column to 
be set into motion with heavy lateral oscillations. 

On the other hand, let us assume the steady load P to be a compression of twice 
the Euler value, while the harmonic part of the load has a frequency f = 2fo€o'/? and 
an amplitude such that h=H/Pz=3.1. We find in this case: 


a, = — 0.25, Bi = 0.775, 
a= 2.00, B. = 3.10, 
a= 15.75, Bs = 6.975, 
a, = n*(n? — 2)/4, B, = 0.775n?. 


We can readily convince ourselves that all points (an, 8,) lie in the stable region of 
Fig. 2. The points (a:, 8:) and (ae, 82) are stable, as one sees from Fig. 2 and the table 
of values of a and 8 for points on the boundary curves given at the end of the paper. 
(Note particularly the values of a and 8 on Cy and C,; for ae —0.25 and the values on 
C2 for a~2.0). The numbers a, can be written in the form a&,=(n?—1)*/4—1/4 
=k?/4—1/4, with k=n?—1; in other words the abscissae a, lie always a distance 
1/4 to the left of the points (k?/4, 0) where the boundary curves delimiting the 
stable regions cross the a-axis. The points (a@,, 0) for m>1 are therefore stable points. 
Also, for 8 not too large the boundary curves lie to the right of the straight lines 
a=k?/4, as one sees from (4.6). Hence all points (a@,, 8,) will be stable if each 8, 
is not too large in comparison with a@,, and this condition is certainly fulfilled in our 
case for n=>2. Note, for example, that 8 must be taken larger than 8 for a point of 
instability when a=8.75 (that is, a value 1/4 less than 9). For ag=15.75 we have 
B; only 6.975 in value so that (a3, 83) is certainly stable. Since the a, increase like m* 
while the 8, increase only like n?, it becomes obvious that all (an, B,) are stable. 
The column is therefore stable even though the steady value of the load is twice that of the 
Euler load.'° However, the total compressive load always, as in this case, drops below 
Pz in value during at least part of the cycle if the column is stable: we have seen 
(cf. 4.5)) that the inequality an+|Bn| >0 holds for stable solutions; in particular, 


for n=1 this leads to 


10 A. Stephenson [13] appears to have been the first to point out the possibility of such phenomena 


in general. This paper appeared in 1908. 
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p+|h|>—1, (5.1) 


as one sees from (3.4) and (3.5), and our statement follows from (2.13). 

Thus there exist both stable and instable sets of values (p, h, f). However, our 
definition of stability leaves out of account a possibility which is always inherent in 
any physical problem, i.e., that slight changes in the parameters of the problem 
(p, h, and f in our case) may be sufficient to cause a stable motion to become an in- 
stable one. A sel of values (p, h, f) should be considered stable in any proper physical 
sense only if a complete neighborhood of these values exists which is made up entirely of 
what we have defined as stable sets of values. 

We proceed to show that the problem of the column never has a stable solution 
in this more restricted sense; i.e., we show that arbitrarily small changes 6f in f and 
5p in p, for example, can always be found such that (p+), h, f+6f) is instable no 
matter what values are chosen for p, h, and f. This is done by showing that a certain 
pair of values (a@,, 8.) becomes instable when properly chosen but arbitrarily small 
changes are made in f and p. Our statement follows from (3.4) and (3.5) and the 
character of the instable regions of the Mathieu equation for high values of a. We 
write equation (3.4) in the form 

atl*/n = (eofe/P)"*(1 + p/n?)s, (S.2) 
and show first that this equation can always be satisfied by taking for a, the square 
of an integer, provided only that f is changed by a small amount 6f and 1 is a suffi- 
ciently large integer: the real number (eof, /f?)"4 can be approximated as accurately 
as desired by a rational number N/n. It is clear that m can always be chosen so large 
that an arbitrarily small change 6f in f will suffice to make the right hand side of 
(5.2) exactly equal to N/n. Hence a,=N‘ and our statement is proved. It is also 
evident that an a, of the form n?/4 could have been determined in the same manner. 
We have thus determined a point (a, 8) for which a=n?/4, n and f+6f being now 
considered as fixed. We recall the fact that the instable regions of the Mathieu 
equation cross the a-axis at right angles at the points where a=n?/4 and that these 
regions for high values of m are narrow strips which remain (for not too large values 
of 8) very near to the vertical straight lines a=mn?/4. Since the values of 8, increase 
like n®, while those of a, increase like n‘ it becomes evident that a small change 6p 
in the value of ~ in (3.4) will be sufficient to cause the point (a@’, 8) corresponding to 
the values p+4, h, f+6f to fall inside an instable region of the Mathieu equation. 
We repeat: no values of p, h, and f (h, f#0) can be found such that the column is 
stable when small variations in these quantities are permitted. 

In the actual physical problem, however, there is an important element present, 
i.e., viscous damping, which has been neglected so far. In a later section we shall 
show that the presence of even the slightest amount of viscous damping will suffice 
to make all values (a, 8) stable for which a2 ay>0, and | B| <a, when aq is a certain 
constant which may be large. In other words, damping acts in such a way as to cut 
out the narrow instable strips which occur for large @ in the regions for which 
| B| <a. Under these circumstances it becomes sufficient to test only a certain finite 
number of the points (a@,, 8,) for stability. Thus the column may be stable if viscous 
damping is present even when small variations in the quantities p, h, and f take place, 
though, as we have seen, this is not the case without damping. 
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Figures 3, 4, and 5 show the stable values of f and h/ (frequency and relative 
amplitude h=H/Pz of the vibratory part of the load) for the values p=P/Pr 
= —1.5, —1.0, and 1.0 respectively. The stable regions are shaded.'! These diagrams 
have been constructed on the assumption that the amount of viscous damping is 
large enough that values of.@ larger than 10 can be ignored. In other words, Figs. 
3, 4 and 5 were constructed by combining the stability regions of Fig. 2, which in- 
cludes values of a up to 10 only, for a suitable number of values of n. 











Fic. 3. 


The general character of Figs. 3 and 4 is typical for the cases in which p< —1, 
i.e., in which the steady part of the load is a compression larger than the Euler 














load. We note that the shaded stable regions for p= —1.5 are much smaller than 
p=-! 
st 
h 
2+ 
i- 
0 2 4 ~ 8 10 2 
2 
fob, 
Fic. 4, 
those for p= —1.0, as was to be expected; for the higher values of the steady 


compressive load beyond the Euler load it is necessary to make more accurate 


11 Without damping, as we have seen, there could be no stable regions though there are stable points. 
It would have a certain mathematical interest to investigate the set of stable points in detail in this case. 
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adjustments in the frequency and amplitude of the oscillatory part of the load in 
order to obtain stability. The full lines which cut through the shaded regions in the 
figures are not really curves; they represent, rather, narrow instable regions. How- 
ever, the two curves in Figs. 3 and 4 which appear to be straight lines running 
near the h-axis indicate narrow stable regions. Fig. 5 is typical for all cases in which 
p> -—1.0, i.e., for cases in which the steady part of the load is either,a tension or 
a compression less in value than the Euler load. In these cases the column is stable 
for all frequencies when h =0; it is in fact stable almost everywhere in the neighbor- 
hood of the axis h=0. 
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It is of some interest to consider the special case in which the amplitude H of the 
oscillatory part of the applied load is very small so that the values of 8, are small 
(for m not too large). We note that the natural frequencies f, of the free lateral oscilla- 
tions of the rod under steady load (that is, in this case, for H=0) are given by 
fn=fa,"'/? as one can readily verify. From Fig. 2 we observe that the rod is instable 
for small values of 8 when a, =k*/4, k being any integer. Hence instability occurs for 
small amplitudes of the oscillatory part of the load whenever 


f = 2fn/k, k= P 2, 3,°°° ’ (5.3) 


that is, whenever the load frequency is twice any integral submultiple of a natural 
frequency of oscillation. At such frequencies one could expect that heavy oscillations 
would be built up."* However, the most favorable case for the production of oscilla- 
tions is, in general, that for which n»=k=1. Consider, for example, the case p=1. 
For n=k=1 we find readily that f/foeo!/?=8?=2.83, and one readily sees from 
Fig. 5 that this furnishes the most favorable frequency for instability at small am- 
plitudes of the oscillatory force. 

6. The flexible string under harmonically varying tension. With only slight modi- 
fications our preceding results can be used to discuss the problem of the vibrating 


2 This problem has been considered both experimentally and theoretically by I. Utida and K. Sezawa 


[16]. 
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string subjected to a harmonically varying tension.’* We have only to set J=0 in 
(2.1) to obtain the fundamental differential equation. The tension F(¢) in the string 
is assumed given by (2.2) and the same independent variables as before are intro- 
duced. However, the parameters p and h in (2.3) can obviously not be used here. 
Instead, we introduce the quantities 


b = P/EA, kh, = H/EA. (6.1) 


We may assume for w the expansion (3.2) for a string with fixed ends and will 
obtain (3.3) as differential equation for the quantities F,(%) if we now define a, and 
8, by the equations 


a, = n*f*p./f, 8, = n°frh,/f. (6.2) 


The investigation of stability involves the same ccns‘derations as for the column, 
and much the same general remarks might be made as were made in the case of the 
column. For example, if P>0 and P>|H), i.e., if the force applied to the string is 
never a compression, viscous damping acts in such a way as to cut out the instable 
regions of Fig. 2 for sufficiently large values of a. Hence it is possible to construct 
a diagram for the determination of the stable values of p, h, and f in the same manner 
as for the column. Figure 6 shows the stable regions (shaded); the quantity f/fop,)/? 
is taken as abscissa and H/P =8,/a,=h,/p, as ordinate. 
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It is readily seen that the natural frequencies f, for the free lateral oscillation of 
the string (under constant tension) are given by f, =fa,'/, just as in the case of the 
column. The string is instable for low amplitudes of the oscillatory part of the tension 
when a, =k?/4, k=1, 2, 3, +--+. In this case we know in addition that f, =f, in 
which f; is the fundamental frequency of the string. Hence “resonance,” that is, 
heavy oscillations for low amplitudes of the applied oscillatory force, will occur when 





% This problem was first discussed by Lord Rayleigh [11]. The problem was discussed later by 
A. Stephenson [12], and [13], and by C. V. Raman [10]. 
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f= 2nfi/k, mn, k =1,2,3,-°-, (6.3) 


that is, at twice any rational multiple of the fundamental frequency of free lateral 
oscillation of the string. However, the most favorable case for the production of 
oscillations is readily seen to be that for which n=k=1 (i.e., that corresponding to 
f/fop."!*=2.0). In Melde’s experiment lateral oscillations of a string are produced 
in accordance with (6.3) by attaching one end of the string to the prong of a tuning 
fork. 

There is one marked (though not unexpected) difference between the behavior 
of the column and that of the string: it could be shown that the string is never stable 
even with viscous damping if the load on it becomes a compression during any part 
of the cycle. For stability of the string we must always require P2 | | : 

7. The effect of damping. If it is assumed that there is a lateral damping force 
acting on the column that is proportional to the velocity dw/dt, the differential equa- 
tion (2.1) is readily seen to be modified by the addition of a term 6(0w/dt), 6>0, to 
its left hand side. With the same notation as before we find as differential equation 
for the functions F,(@): 

nN Mn oF =0 7.1 
3% r= + ethene , (7.1) 
where 


= 5/4rm/f, (7.2) 


and subscripts have been dropped. 

The general theory of equation (7.1) could be developed in the same way as that 
for the Mathieu equation without damping (for a treatment which includes a damp- 
ing term, see the papers of G. Gorelik [3]). In particular, the a, B-plane could be 
divided into stable and instable regions. We confine ourselves here to one special 
problem, i.e., to a discussion of the behavior of the solutions of (7.1) for a given 
value of v and large positive values of a. We assume also that |8| <a. 

Upon making the substitutions 


F = «~"G, a’ =a-—y (7.3) 
Eq. (7.1) becomes 

. (a’ + B cos 9)G = (7.4) 

di? 


Obviously, if G is bounded, F is not only bounded but approaches zero as # increases. 
Also, even at boundary points (a’, 8) separating stable and instable regions of 
(7.4), the corresponding solutions F tend to zero since no solution G of (7.4) in- 
creases faster than # in this case. If the amount of damping is slight (that is, if » 
is small), the boundary curves for (7.1) would lie near those for y=0, but they would 
not intersect the a-axis except at the origin since all solutions F of (7.1) are clearly 
bounded for 8=0, a>0. This reasoning makes it seem rather evident that the nar- 
row instable regions which occur for large positive values of a when |8| <a are cut 
out when a damping term is added. 

We proceed to give a proof of the following statement: if v>0 and 'a| <a, all 
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solutions of (7.1) are stable for all values of a which exceed a certain value ao>O. It 
was pointed out earlier that there exist two linearly independent solutions G; and G2 
of (7.4) such that Gi@+27)=KGi@) and KG:@0+27) =G.(@) with | K|>1 in case 
(a’, 8) is in an instable region for (7.4). We know also that K is a real number in this 
case. The solutions of (7.1) will, however, remain bounded even in such an instable 
case for (7.4) provided that 

ery >| KI, (7.5) 


as one sees from (7.3). Consequently our statement will be proved if we can show 
that | K|—1 as a’—~o. This we prove through the use of the following asymptotic 
formula for K, valid under our assumptions, which has been given by Strutt [14]: 


K + 1/K = 2 cosh X cos ¢ + O(1/V/a’) (7.6) 


in which 
t+rA/-1= f (a’ + B cos 8)dd, (7.7) 
0 
and O(1/\/a’) means that all terms neglected are of order 1/\/a’ or higher. Since 
we assume that |8| <a’ the integral in (7.7) is real and \=0. We have, therefore: 


|K +1/K| <|2cos¢| + O(1/Va’) < 2+ O(1/Va’). (7.8) 


Since K is real it is readily seen that 
2<|K+41/K|, (7.9) 


equality holding only for | K| =1. From this and inequality (7.8) it follows at once 
that 
|K|—1 when a’>o, (7.10) 
In the case of the column we note from Eqs. (3.4) and (3.5) that |B..| <a, for 
sufficiently large m and that a,—« with m. The assumptions under which (7.10) 
was derived are thus fulfilled in this case. When damping is present we are there‘ore 
justified in neglecting all values of a larger than a certain positive value ap in discuss- 
ing the stable values for the column. Our diagrams were drawn under the assumption 
that a)»=10. In the case of the string, a, and 6, increase at the same rate with in- 
crease of 2; consequently our conclusions regarding the effect of damping in this 
case are valid only when P>|H| (which ensures that | Bn| <@Q@,,) i.e., when P isa 
tension and H is such that the total force in the string is always a tension. 


APPENDIX 
Sufficient conditions for stability. For stability we required always that the solu- 
tion 
w= > F,(d) sin nt (A1) 


n=1 


of our problems be bounded for arbitrary initial conditions; it is thus necessary to 
assume for stability that each F,@) be bounded for 0S3<@ (@ is essentially the 
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time variable). In this appendix we prove a statement made at the end of section (3) 
to the effect that the series will converge for all # if it converges for d =0 and if the 
F,,@) are all stable solutions of the Mathieu equation. 

In order to state our theorem precisely we introduce the series 


} D,, sin né, > V, sin né (A2) 
in which D, and V, are defined by 
dF, 
D, = F,(0), V2. = —| . (A3) 
dd \s=-0 


We assume that the series (A2) are such that 
Dd {rn 


in which p, is a certain positive quantity and a, is one of the two parameters in the 
Mathieu equation for the functions F,(@): 





+ Patin | Vil} < @, (A4) 





D, 





1F,, 
—— + (a, + B, cos ¥)F, = 0. (AS) 
dd? 


We assume in addition that the F,,@) are stable solutions of (AS) for which 
|Bn| < kan, OS k<1, (A6) 
at least for all »>N, say.“ Under these assumptions we show that: the series 


dF,(8) . 
— sin n& 
v 





bm F,,(8) sin nE and ys 


n n 


converge for 0<3 < « in the same sense as the series (A2), 1.e., the convergence relation 





_y2| dF, (0 
> {en| Fad) | + prota | } < © (A7) 


holds for OS 8< ~, 

If it were assumed that p,=1 in (A4) then > .F.@) sin n— would converge, but 
its derivative with respect to? would not necessarily converge. If p, were assumed 
to be al”, the differentiated series would converge. In our cases al’? is of order n 
for the string and of order mn? for the column. To assume p,=a,” in (A4) would 
therefore not seem unduly restrictive when it is considered that the series (A1) 
should be assumed to converge when it is differentiated twice with respect to & in 
the case of the string and four times with respect to in the case of the column. 

We prove our theorem by showing that every stable solution of the Mathieu 


equation 


d°*F 
Q + (a+ Bcosd)F = 0 (A8) 





14 These latter conditions are fulfilled in the stable cases for both column and string. This follows from 
(3.4) and (3.5) for the column, and from (6.2) and the fact that | he| <p, in the case of the string. 
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for which 
\B|<ka, OSk<1, (A9) 
and 
F(0) = D = V (A10) 
*( = ; coientaaiaiiaas = f 
dd \0=0 ' 
satisfies the inequality 
| dF (8) | | 
| Fo) | + ~~ <C{|D| +a7]|V|}, (A11) 
| ¢ | 


for 0<0 < «, C being a constant which depends only upon k. Upon reintroduction 
of the subscript in (A11) followed by multiplication with p,>0 and a summation 
with respect to , it is clear that (A7) would result from (A4), since C is independent 


of x. 
We proceed to establish the validity of the inequality (A11). For this purpose it 


is convenient to introduce a new independent variable ¢ as well as a new dependent 


variable f in (A8) as follows: 


3 
o(0) = f x)/*dd, x =a-+tfBcos#, (A12) 
0 


f= yl/4P, (A13) 
J 


In these variables the differential equation (A8) becomes 


d*f Bcosd’ - 56? sin? 0\ | 
— -+ (1 + — + — f= 9, (A14) 
4x? 16x? 








dg? 


or, as we prefer to write it 








af 
ier" (A15) 
dg? 
with 
Bcosd? 568? sin? # 
y=- a(= tod “), (A16) 
4x? 16x? 


From now on we consider f(g) to be the solution of (A15) which satisfies the 


initial conditions 
df aes 
—=1 l — 1. (A17) 


f(0) = 1, 1, ~=V\ 
dy 


It is then readily verified that f(¢) and its derivative satisfy the integral equations 


i id ° 
f(y) = ef? —— (< f y(r)f(r)e—*dr — ci f vor)f)edr), (A18) 
2a 0 0 


df(¢) : a earner ’ 
—_——— ie‘? + — (cr f y(r)f(r)e~*dr — e ef virifr)erdr). (A19) 
0 0 


dy 2a 


15 This transformation is frequently used in the treatment of various questions relating to the asymp- 
totic behavior of the solutions of certain types of second order ordinary differential equations. 
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From the general theory of the Mathieu equation it is known that every stable 
solution F@) of (A8) can be expressed in the form H@)e, in which H@) is a periodic 
function of period 27 and a is a real constant. It follows from (A13) that f(¢@)) 
can be expressed in the form h@)e“ with h = Hx!/*; h@) is thus also periodic of period 
2m ind’. Consequently we may write 


G = max | f(8) | = max | (8) | = maxjjse| 2) | = maxjsjse|f()|. (A20) 


The validity of (A20) is the essential point in our proof; because of it, bounds for our 
quantities in the interval —7 <0 Sr hold also for OS 0< &. 
We find from (A16), the definition of x in (A12), and (A9) that 


ly |S &/4(1 — B*) + 5k2/16(1 — A) = 7. (A21) 
We note also that 
o(r) SrVat+BsrVaVithk, (A22) 


as one sees from (A12). Finally we obtain from (A18) the following inequality for 
G=max |f(@)|: 
IinrV1 +k 
G s 1+——__G. (A23) 
Va 
In view of our purpose it is permissible to assume from now on that 
a => ap > F*x*(1 + 2k) = a; (A24) 


once this is done (A23) may be written in the form 


G < 1/(1 — Van/ao) = Go. (A25) 
In a similar fashion we can show that 
d 
max Ae) < Go, (A26) 
dy 








since df/dg satisfies (A19) and, like f(g) itself, can be written in the form hoe 


with h of period 27 inv. 

Since the function f(¢@)) given by (A18) and its complex conjugate are linearly 
independent solutions of (A15) it follows that we may write the general real solution 
F of (A8) in the form 


F(8) = Re Cx~—"/4f(8), C=A — iB, (A27) 
in which Re means that the real part of what follows is to be taken, and A and B 


are real but otherwise arbitrary constants. The quantity dF/dd is then given by 


the expression 





dF f _— | (A — iB)f(e(8))6 sin “ 
a SS 2 a ys A- —_ . A28 
Re \( + B cos 8)(A — iB) rm + 1/4 (a 4 B cos 0)e (A28) 


We find at once, since 
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1 
g(0) = 0 and sf =1, 
dy g=0 
dF 
D =F(0) = (a+ 6)""44, Ve= sacl = (a + B)'/4B, (A29) 
dd \o~0 
from which we obtain 
A=(a+ ,)'4D, B=(a+ 8)”. (A30) 


For | F@)| we then have the inequality 


1+ k\'/* 1 
\(—) | D| + —_—__| v | a 
1 -— oet/2(1 — 2) 1/4 


V | (A31) 


’ 
| 





F(@) | 


IA 


< Po| D| + goa*/? | 


in which fp and go depend only upon the constant k introduced in (A9), and Gp is 
the bound for max | f(@)| given in (A25). From (A28) we find 


dF | k aie ae 
— | < 1/4 —— (a — B)-"/4| C|Go + (a2 + B)"*| C| Go 
| dd | i-—k 
< a'/p,|/D| +q| VI, (A32) 
where 
jl k (- +" : a | 

= ff ————} + (1 + 9G (A33) 

Ps 4a? 1— kM — i a 


and q; is of similar nature. The quantities p; and qi, like po and go in (A31), depend 
only upon k. Division of both sides of (A32) by Va, followed by addition to (A31) 


yields 
| dF | 
|F| +a? ry < p|D\| +a—g| VI, (A34) 


‘eG. 


which establishes the validity of (A11) and thus completes the proof of our theorem. 


Coordinates of Points on the Boundary Curves of Fig. 2. 





B a(Co) a(C,) a(.S;) a(S) 

0.0 0.00000 0.25000 0.25000 1.00000 
0.2 —0 .01966 0.14525 0.34475 0.99667 
0.4 —0.07510 0.03191 0.42796 0.98670 
0.6 —0.15836 —0.08872 0.49816 0.97018 
0.8 —0.26148 —0.21555 0.55906 0.94724 
1.0 —0.37849 —0.34767 0.59480 0.91806 
1.2 —0.50535 |  —0.48430 | 0.62006 0.88284 
1.4 —0.63942 —0.62480 | 0.63015 0.84183 
1.6 —0.77898 |  —0.76867 0.62592 0.79529 
1.8 —0.92281 —0.91545 | 0.60857 0.74349 
2.0 —1.07013 —1.06480 0.57950 0.68672 
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a(Co) a(C,) a(S1) a(S2) 
—1.22031 —1.21640 0.54012 0.62526 
~1.37291 —1.37002 0.49174 0.55938 
—1.52760 —1.52544 0.43554 0.48935 
—1.68410 —1.68248 0.37253 0.41542 
—1.84221 —1.84098 0.30357 0.33785 
—2.00175 —2.00081 0.22938 0.25684 
—2.16258 —2.16185 0.15057 0.17263 
—2.32457 —2.32402 0.06763 0.08541 
—2.48764 —2.48720 —0.01901 —0.00468 
—2.65168 —2.65134 —0.10899 —0.09734 
—2.98242 —2.98220 —0.29781 —0.29009 
—3.31627 —3.31614 —0.49688 —0.49171 
—3.65286 —3.65277 —0.70474 —0.70124 
—3.99186 —3.99180 —0.92026 —0.91787 
—4.33302 —4.33298 —1.14253 —1.14088 
—4.67611 —4.67609 —1.37085 —1.36970 
—5.02097 —5.02096 —1.60460 —1.60383 
—5.36744 —5.36743 —1.84328 —1.84271 
—5.71537 mam 4.) . —2.08644 —2.08607 
—6.06467 —6 .06466 —2.33382 —2.33353 
—6.41522 —6.41522 —2.58498 —2.58478 
—6.76694 | —6°76694 —2.83970 —2.83955 
—7.11974 —7.11974 —3.09772 —3.09761 
—7.47357 —7.47357 —3.35883 —3.35875 
—7.82835 | —7.82835 —3.62283 —3.62277 
—8.71911 | —8.71911 —4,29436 —4.29434 
-9.61474 | —9.61474 —4.98065 —4.98064 
—10.51465 | —10.51465 —5.67983 —5.67982 
—11.41834 —11.41834 —6.39044 —6.39043 
—12.32542 —12.32542 —7.11126 —7.11126 
a 

—13.23556 | —13.23556 —7.84129 —7.84129 
—15.06389 —15.06389 | 9.32566 —9 32566 
~16.90154 | —16.90154 ~10.83807 —10.83807 

a(¢ 2) a(C3) a(S3) | a(S,) 
1.00000 2.25000 2.25000 4.00000 
1.01633 | 2.25225 2.25275 4.00133 
1.06171 | 2.25808 2.26203 4.00530 
1.12806 2.26622 2.27933 4.01181 
1.20733 | 2.27554 2.30589 4.02075 
1.29317 | 2.28515 2.34258 4.03192 
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8 ae(Cs) | a(C3) | ae(Ss) a( Si) 
1.2 1.38126 2.29429 2.38967 4.04512 
1.4 1.46860 2.30233 2.44680 4.06010 
1.6 1.55305 2.30878 2.51308 4.07660 
1.8 1.63302 2.31323 2.58723 4.09433 
2.0 1.70727 2.31536 2.66777 4.11301 
2.2 1.77487 2.31495 2.75314 4.13236 
2.4 1.83509 2.31175 2.84194 4.15212 
2.6 1.88745 2.30568 2.93284 4.17199 
2.8 1.93163 2.29660 3.02467 4.19175 
3.0 1.96752 2.28448 3.11640 4.21115 
3.2 1.99517 2.26925 3.20712 4.22997 
3.4 2.01478 2.25092 3.29604 4.24800 
3.6 2.02665 2.22950 3.38247 4.26507 
3.8 | 2.03118 2.20500 3.46578 4.28099 
4.0 2.02881 2.17748 3.54547 4.29563 
4.4 | 2.00521 2.11356 3.69216 4.32053 
4.8 1.95947 2.03826 3.81969 4.33886 
5.2 | 1.89487 | 1.95216 3.92636 4.34996 
5.6 1.81419 1.85589 4.01149 4.35338 
6.0 1.71968 1.75014 4.07538 4.34881 
8.0 1.09281 1.09947 4.12172 4.20467 

10.0 0.28857 | 0.29018 3.84895 3.87349 
12.0 | 0.63494 | 0.63452 3.38071 3.38817 
14.0 | —1.64702 — 1.64690 2.77777 2.78016 
16.0 | —2.72859 |  —2.72855 2.07287 2.07367 
18.0 |  —3.86669 | 3.86668 1.28641 1.28668 
20.0 | —5.05198 | 5.05198 0.43241 0.43251 

B | (C4) a(Cs) a(Ss) a(Ss) 
0.0 4.00000 6.25000 | 6.25000 9.00000 
0.2 4.00134 6.25083 | 6.25083 9.00057 
0.4 | 4.00538 6.25333 6.25333 9.00229 
0.6 | 4.01226 6.25750 6.25751 9.00515 
0.8 4.02215 6.26334 6.26337 9.00915 
1.0 | 4.03530 6.27084 6.27094 9.01430 
1.2 4.05204 6.27999 6.28025 9.02060 
1.4 | 4.07273 6.29077 6.29134 9.02806 
1.6 4.09776 6.30317 6.30427 9.03667 
1.8 | 4.12755 6.31714 6.31911 9.04643 
2.0 | 4.16245 6.33264 6.33594 9.05735 
2.2 | 4.20283 6.34961 6.35487 9.06943 
2.4 4.24889 6.36800 6.37604 9.08267 
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8 | a(C.) (Cs) a(Ss) (Se) 
2.6 4.30085 6.38773 6.39956 9.09705 
2.8 | 4.35867 6.40871 6.42560 9.11259 
3.0 | 4.42220 6.43085 6.45432 9.12927 
3.2 4.49121 6.45406 6.48591 9.14707 
3.4 4.56533 6.47821 6.52052 9.16600 
3.6 | 4.64406 6.50321 6.55837 9.18603 
3.8 | 4.72688 6.52893 6.59962 9.20714 
4.0 | 4.81318 6.55525 6.64444 9.22930 

| 
4.4 4.99383 6.60921 6.74533 9.27671 
4.8 5.18127 6.66411 6.86185 9.32798 
5.2 5.37113 6.71898 6.99394 9.38281 
5.6 5.55951 6.77289 7.14093 9.44078 
6.0 | 5.74803 6.82500 7.30201 9.50150 
8.0 6.50217 7.03409 8.23272 9.82875 
10.0 6.89864 7.11706 9.16125 10.14742 
12.0 6.97136 7.05384 9.87814 10.40143 
14.0 6.82083 6.85144 10.30874 10.55621 
16.0 6.51561 6.52721 10.48838 10.59848 
18.0 6.09463 6.09902 10 .48167 10.52959 
20.0 | 5.58132 | 5.58302 10.33749 10.35813 
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ON MOMENT BALANCING IN STRUCTURAL DYNAMICS* 


BY 
R. E. GASKELL** 


Brown University 


1. The method of moment balancing. In recent years several writers in this coun- 
try have developed the method of moment balancing in the analysis of continuous 
beams and frameworks. Mention should be made especially of the basic paper by 
Hardy Cross.'? One could also classify as related procedures the method of balancing 
angle changes given in a paper by L. E. Grinter,* and the whole field of relaxation 
methods being investigated by R. V. Southwell.‘ That such interest is taken in these 
methods would seem to indicate that their extension to the dynamics of beams and 
frameworks might be desirable, and it is the purpose of this article to provide at least 
the beginning of this extension. 

We assume that we are dealing with plane structures on which loads are acting in 
the plane of the structure. Members of the structure consist of uniform straight 
beams; and they meet in stiff joints, which are assumed to be fixed against translation. 
All connections to a foundation are either built-in or hinged. 

The method of moment balancing depends upon three very simple ideas, namely, 
fixed-end moment, stiffness and carry-over factor. We give their definitions here: 

The “fixed-end moment” at the end of a member is the moment which would 
exist at that end if all joints to which it is connected were fixed against rotation. 

If one end of a member is simply-supported, its “stiffness” is the moment re- 
quired to produce unit rotation of that end. The other end may be built-in, simply- 
supported or free. 

The “carry-over factor” is the numerical value of the moment induced at one 
end of a member by a unit moment acting at the other end. 


Methods of finding these characteristics of beams and other components of a 
structure are numerous and well-known. Having determined them for all components 


* Received Feb. 3, 1943. 

** This paper was prepared under the direction of Professor W. Prager, whose helpful suggestions 
and valuable assistance are gratefully acknowledged. The author is a Fellow under the Program of Ad- 
vanced Instruction and Research in Mechanics at Brown University. 
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of a framework, we assume that all joints of the framework (in Fig. 1, for example) 
are fixed against rotation; and determine the resulting fixed-end moments acting at 

the ends of each member. Built-in, simply- 
SS ¢ supported or free ends are not considered as 
joints. Then at any joint, say D, a moment 


A D Pad - Ae equal but opposite in sign to the sum of its 
+ za 





fixed-end moments is applied, representing the 
effect of releasing the joint. This moment is 
RSS 2 distributed to the members AD, CD, BD, ED, 
meeting at D, in proportion to their stiffnesses, 

Fic. 1. since all members meeting at D rotate through 

the same angle. The share falling to each mem- 

ber is called the “balancing moment” acting at the end D of this member. The joint 
D is now balanced, but the new balancing moment Mp, acting at the end D of AD 








A8 


will induce an additional moment 


Map = CapMpa 


at the opposite end A. Cy p is the carry-over factor for the member AD, and the mo- 
ment Mp, is said to be “carried over.” Likewise, moments are carried over to C 
and E, but none to B since Cpzg=0. The joint D 

is again locked—this time in its balanced posi- B = C 


tion—and the process repeated for all joints of 








the framework until the balancing moments are d , a 
negligible. The order of choosing unbalanced 
joints for balancing is not obligatory, but usu- 
ally the joint with the largest total unbalanced 

moment at any given stage is balanced. Signs / 
of the moments are chosen so that a positive 
moment acting on the end of the beam tends to 











rotate it in a clockwise direction. Likewise, a 
rotation in the clockwise direction is consid- LL wey? 
ered positive. 

Example 1. As a simple example consider the Fic. 2. 


rectangular bent formed of uniform and equal 
bars, illustrated in Fig. 2. All of the bars are of equal stiffness and the carry-over fac- 
tor in each case is 1/2. The only non-vanishing fixed-end moments are —.125Pl 
and .125P/ at the left and right ends of the horizontal bar. The calculations used 
in the method of moment balancing are shown in Table I. In a given column, say 
that headed Mcz/Pl, we find recorded successively the fixed-end moment and the 
balancing moment. These are added, and since at this stage Mcz+Mcp=0, the joint 
C is balanced. The balancing moment has been carried over to column Mgc/Pl, and 
the joint B is balanced next. The same steps are followed until after five balancings 
the moments to be carried over are negligible. The results obtained agree with those 
computed by other methods. 

2. Dynamics of a simple beam. It is clear that if we can set up analogous defini- 
tions for fixed-end moment, carry-over factor and stiffness for a beam on which an 
oscillating force is acting, and if we can find these characteristics for the oscillating 
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beam, it may be possible to use the method of balancing moments just as it is in the 
dynamic case. A procedure adapted to this purpose can be found in an article of 
W. Prager’s,® the essentials of which will be given here. 







































































TABLE I. 
B Cc 
Mas/Pl Mepa/Pl Msc/Pl Mcs/Pl Mep/Pl Mpc/Pl 
.000 .000 —.125 .125 .000 .000 
— .063 — .062 
— .032 — .031 
.000 .000 —.157 .062 — .062 —.031 
.078 .079 
.039 .040 
.039 | .078 — .078 .102 — .062 — .031 
— .020 — .020 
— .010 — .010 
.039 .078 — .088 .082 — .082 — .041 
| .005 .005 
.002 .002 
.041 .083 — .083 .084 — .082 — .041 
— .001 — .001 
.041 .083 — .083 .083 — .083 — .041 


























The differential equation for the deflection, y(x, ¢), of a uniform beam with no 

external load is taken as 

d°y dty 

p—— + EI — = 0, 

ot? Ox‘ 
where yu is the mass per unit length of the beam and EJ 
is its flexural rigidity (Fig. 3). Following a well-known @ Ely d 
procedure we write y(x, t)=u(x) cos wt, u(x) being the _ 
amplitude of the assumed harmonic motion and w its ¥ 
circular frequency. Hence —" 

d‘u + a! 9 

—-+n‘u = 0, 

dx* 





where n‘=w*u/EI; and from this equation 
u(x) = A cosh nx + B sinh nx + C cos nx + D sin nx. 
It is convenient to express the four constants of integration in terms of four quanti- 


ties of immediate physical importance: the amplitudes of the moments acting on the 
ends of the beam, and of the displacements at the ends of the beam. This can be done 


by use of the relations 


’ W. Prager, Die Beanspruchung von Tragwerken durch schwingende Lasten, Ingenieur-Archiv 1, 527 
(1930). 
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uy = A + GC 
u, = AcoshA + BsinhrA + C cosdA + Dsin\, 
_ au 
M)=—EI—| =-—A-—C}EIn},{ 
dx? z=0 
d*u : 
M, = EI i = {A cosh\ + B sinh A — C cosd — D sind} Eln’, 
AX" 5 cal 


where \=m/. It can be seen that the amplitudes of the deflection, angle of rotation, 
bending moment and shear for any value of x will involve linearly the amplitudes of 
the end deflections and end moments. These quantities can be expressed in much 
simpler form if the following functions and abbreviations are introduced: 


ae r 
o(A) = (coth A — cot A)/2\, o(A) = — (coth A + cot dA), 


ll 


< r 
¥(A) = — (csch A — csc A)/2A, ¥(A) = — (csch A + csc A), 


VY = W/E. 
Then we find the following expressions for the amplitudes of the angles of rotation at 
the end points (Fig. 4): 





Ftc. 4. 
f uob(r) u(r) 
us = ——— + —— + Mdl'6(0) — MV), (1) 
uow(d) uip(r) 
fm: - -— Mil'¥(d) + Mil’9(); (2) 


and for the amplitudes of the reactions: 


4 


M M, 


arr r 4 Q0O— a 
Ro = — Eluy = uUp— G(r) + uw, —VA) -- (A) - —yvA), 3 
' 0a 90) + 1 VQ) — = 80) — =F) (3) 
— ae ae Mo - Mi_ 
R, = = EIu, = — uo ry (A) — @e PY (A) samed 7 ow — “yo. (4) 


If the beam, simply supported at both ends, is loaded at its center by an oscillating 
load, P cos wt (Fig. 5), 
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l l 
My = % = wv’ (5) = 0, R (5) = P/2. 
Z 2 


Then, from (2) and (4), 





2 PX al & I’ l ry Pp 
"(x )8(z)+5 ¥G)() = J | d 
: ae a 
ik tae L, 2 (4)s(2 P j—l 
— 7 alae Male Mapeer ae ne 
so that 
l 
u(> = PI*l/ (a), (5) 
l ‘ 
M (=) = — Pl&(a), (6) 
where 
@(\) = — (tanh $d — tan $d)/4A, #(\) = (tanh 3 + tan 3A)/4A. 


Also, from formulas (5) and (6), and with (1)—(4), we find that 


ug = PI'¥(d) = — ui, (7) 
and Se 
Ro = P¥(A) = — Ri, (8) 
where 
W(A) = — (sech 3 — sec 3A)4A2, WA) = (sech $+ sec $d)/4. 


Now, if the beam in question is on unyielding supports but has moments acting 
on its ends in addition to the load acting at its center, we find by addition of (1) and 
(7) that 

us = Mol’o(d) — Mil’) + PI'VA); (9) 


and, similarly 


ui = — Mol'W(r) + Mil'o(A) — PU’). (10) 


Obviously, with the equations derived, problems in dynamics of frameworks are 
reduced to problems in statics of frameworks. To facilitate this work, tables of the 
functions (A), dA), WA), JA), PA), BA), WA) and WA) are available.‘ 

3. Dynamic moment balancing. By substituting “moment-amplitude” and “rota- 
tion-amplitude” for “moment” and “rotation,” respectively, wherever they occur in 
the definitions of fixed-end moment, stiffness and carry-over factor, we arrive at 
suitable definitions for the corresponding dynamic quantities. These three quantities 
will give us a basis for the application of the moment balancing method to problems 
in dynamics of frameworks. 


* K. Hohenemser and W. Prager, Dynamik der Stabwerke, Julius Springer, Berlin, 1933. 
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First let us consider the amplitudes of the moments acting at the ends of a cen- 


mig— [gm 





Fic. 6. 
Uy=/ 
Me 
fa) 
oe, )M 
(b) 
Up ———~ai )™ 
kc) 
Fic. 7. 
— Mdl'W(rs) + 


and from (9) 


Moo() — MwA) 


Since M; is by definition the stiffness, K, 
$(A) 


trally loaded built-in beam (Fig. 6). Equa- 
tions (9) and (10) can be applied, and we 
find that 

Moo(d) — My) = — PIV), 
and 

Mw(r) — Mig(d) = — PIV(A). 
Then, since ¢(A) +W(A) =2(A), the ampli- 
tudes of the moment are given by the 
relations 
PI¥(A) 


a Me ee eee 
2@(d) 


Mo (11) 
These quantities give the amplitudes of 
the fixed-end moments for a beam loaded 
at its center with a load of amplitude P. 

The problem of finding the dynamic 
stiffness is illustrated in Fig. 7. If the far 
end of the beam is built-in (Fig. 7a), we 
find from (2) that 


M \l'o(d) = 7 


Il 


AB(A) 





ry Vt [oa) ]}? — va) ]?} 


where 


- ——» 12 
D(A) (2) 


B(A) = cosh X sin A — cos J sinh d, 
D(A) = cosh A cos A — 1. 


Tables exist for these functions and for the quotient B(A)/D(A).® 
For a beam on two simple supports (Fig. 7b), equation (2) gives 


1 = M\l'¢(), 


so 


K = 1/l'¢(0). 


To find the stiffness of the cantilever beam (Fig. 7c), we find from (2) and (3) that 


uo — 
~ 7m + Mil’¢() = 1, 


4 
Uu —_—_— 
” PY 
Hence, 


A4H(A) 


M, 
o(d) — —") =» 6. 





X= Teal — WOR) 


(13) 
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The carry-over factor needs to be found only for the case illustrated in Fig. 7a, 
since it is zero in the other two cases. 
From equation (1) 


Moo(r) — Mi) = 0, 


so that the carry-over factor is defined by 


C=a—=-—- (14) 


We are now equipped to apply the moment balancing method to problems in 
dynamics of frameworks. 

Example 2. Consider again the bent illustrated in Fig. 2, but now suppose that 
the frequency w has a value such that \ =3.30 for each bar. Then we have fixed-end 
moment-amplitudes of —.169P/ and .169P/ at the left and right ends of the horizontal 
bar, equal carry-over factors of 1.22, and equal stiffnesses for each bar. Table II 
gives the calculations involved in solving this problem. The values obtained from the 
12 balancings are correct to two significant figures. 

4. Dynamic balancing of angle changes. The application of the results of Sec- 
tion 2 to L. E. Grinter’s method of balancing angle changes® is not difficult. In 
balancing a given joint, the members of the framework meeting at the joint are as- 
sumed to be simply-supported and disconnected there. Then rotations are forced by 
means of applied moments until the angular discontinuities between the members are 
negligible. To work with rotations rather than moments we require two more defini- 


tions. 


By “angle-change” will be meant the change in slope produced at the end of 


a member either by loads or by an applied end moment. 
The “angle carry-over factor” is the numerical value of the angle change in- 
duced at one end of a member by a unit angle-change imposed upon the other end. 


The amplitudes of the angle changes, at the ends of a simply supported beam, 
due to a central load of amplitude P are seen from (7) to be 


uo =— uy = PIl'V(A). 


The angle carry-over factor can be found by consideration of a simply supported 

beam, one of whose ends is rotated by means of an applied moment-amplitude (Fig. 

8). From equations (1) and (2), uf =—u¢ W(A)/d(A) so that the angle carry-over 

factor is ¥(A)/(A). ‘ ’ 
Similarly, by consideration of equations (1), Ag | 4 u, 

(2) and (4) we can arrive at an angle carry-over 1am 

factor for a cantilever beam: Fic. 8. 


_ GAMA) + MOAWA) - 
[Za) }? — *[9) }? 


Since (0) = J(0) =1, it is seen that this reduces to unity in the static case. 
Continuity is established between a member and a joint by giving the joint a 
rotation-amplitude K,6;/2K, where 6; is the amplitude of the angle change in the ith 
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TABLE II. 




















B ¢ 
Mas/PI Mza/Pl Msc/PI Mcs/Pl Mcp/Pl Moc/PI 
.000 000 — .169 169 | .000 000 
— .084 — .085 
— .103 | | — .103 
.000 000 — .272 .085 — .085 |} —.102 
136 136 
.166 .166 | 
—— ‘ é vs i a | 
.166 136 | —.136 .251 — .085 |} —.103 
— .083 — .083 | 
—.101 —.101 
.166 136 a 293 168 — .168 — 204 
051 .050 
062 062 
| 
228 | 187 =. | 230 — .168 204 
] — .031 — .031 
— .038 ] | — .038 
HT | 
ia nam ¥ oe 1] bas SS Pat 
.228 | .187 — .225 ] 199 — .199 | 242 
019 .019 | 
023 | 023 
- — 1 we ——— a 
251 206 — .206 || .222 — .199 — .242 
mn 049 —.012 | 
— .014 — O14 
mete eee ee int | aortas 
251 206 — .220 ] 211 on 909 | — .256 
007 007 i] | 
009 009 | 
wee ——— \| a are ke | be o 
260 213 =~ 993 | 220 wo S04 | ~ .256 
| — .005 — 004 | 
— .005 i] — .005 
—_ — a oma im 
.260 213 218 215 —, 215 | — .261 
.002 .003 | 
.003 003 | 
263 215 12s | 218 =~ $8 _ 
—.001 | —.002 
— .002 — .002 
263 215 «939 .217 | ST | oo 8 
001 .001 
263 216 = 316 217 ~ 207 — .263 
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member, K; its stiffness, and the summation extends over all bars meeting at the joint; 
and at the same time the end of the member itself is given a rotation-amplitude 
—|0;—K,0;/2K]|. This is done for each member meeting at the joint, thereby bal- 
ancing that joint; then the assigned rotation-amplitudes are carried over, and the 
balancing process continues. 

After the rotation-amplitudes for all joints have been found with the desired ac- 
curacy, the moment-amplitudes can be found from a combination of (9) and (10): 


ud o(d) + ui (A) PI¥(d) 



































M, = —-—— _ , (15) 
Vf [oa]? — Wa)}?} 26a) 
- uo W(A) + ui P(A) PIV () (16) 
4 1 =—_- rine = = Pre es a 
l'{ [ory]? — [Wa)}?} 2A) 
TABLE III. 
e e eee ee 
B C 
6p4/PIl’ @ac/PIl’ | 6cp/PIl’ Ocp/Pil’ 
.0000 .0625 — .0625 .0000 
0357 — .0268 
— .0178 
.0000 0447 — .0268 — .0268 
.0192 — ,0225 | 
.0128 
0192 0192 — .0140 — .0268 
— .0073 0055 
.0036 
0192 .0228 — .0213 — .0213 
0016 — .0020 
.0010 
0208 .0208 — .0203 — .0213 
— .0006 .0004 
.0003 
.0208 .0211 — .0209 — .0209 
0001 — .0002 
.0001 
0209 | 0209 — .0208 — .0209 
— .0001 
.0209 .0209 — .0209 — .0209 
Map= .0418 Pl Mpa=_ .0836 Pl Mce= .0832 Pl 


Mpc= — .0419 Pl Mac= — .0832 Pl Mcp= — .0836 Pl 
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It is interesting to observe that when \=0 equations (15) and (16) reduce to the 
slope-deflection equations for a centrally loaded beam.’ 

Example 3. As an illustration let us solve Example 1 by the method of balancing 
angle changes. For this method, the stiffness of the horizontal bar will be the moment 
required to produce unit rotation of one end while the other end is simply-supported; 
while the stiffness of a vertical bar requires the other end to be built-in. Hence the 
ratio of the stiffness of the horizontal bar to the stiffness of a vertical bar is 3/4. If 
all joints are assumed to be pin-connected, we have angle changes 0gc =.0625 Pll’ 
and 6¢z = —.0625Pil’ due to the load P. For simplicity, u/ is replaced by 6. 

In Table III we find the computation used in solving this example. Joint C is 
balanced first by rotating the member CB through the angle —(1—#)(—.0625) Pll’ 
=.0357Pil’ and the other members of the joint (that is, CD) through the angle 
#(—.0625) Pil’ = —.0268Pil’. Continuity at that joint is then established, but the 
rotation of CD induces a rotation at the other end of the beam, 0g¢ = —.178Pll’. 
This leaves a total unbalance of .0447P/l’ at joint B, which is balanced next. These 
balancings continue until the angle changes to be carried over are negligible. The 
resulting moments, computed from (15) and (16) are also listed, and compare favor- 
ably with the results obtained for Example 1 by moment balancing. 

Example 4. If, now, w has a value such that \=3.00, we find angle changes 
O0ec =.381Pil’ = —Ocz due to the load P coswt. Furthermore the angle carry-over 
factor for the horizontal bar is —.872, and as to stiffnesses, Kpc =.549 1’, Kan =Kep 
= 3.102 1’. Table IV gives the computation involved in 12 balancings of angle changes 
in this case. The values of the moment-amplitudes obtained are compared with those 
obtained by moment balancing. 

5. Convergence of the moment balancing process. Convergence of the process of 
moment balancing can be assured if the frequency of the forced vibration is smaller 
than the first natural frequency of the structure. The first step of the method of mo- 
ment balancing leads to the determination of the amplitudes of the unbalanced 
moments. For the following steps these unbalanced moments are considered as ex- 
terior couples acting on the joints of the structure. In the type of structure considered 
here (joints fixed against translation) the amplitudes of displacement and bending 
moment of any member are completely determined by the frequency w and the 
rotation-amplitudes at the two ends of the member. If a set of values of the rotation- 
amplitudes at the m joints of the structure is assumed, it is therefore possible to com- 
pute the amplitudes of the periodic couples which must be applied to the joints in 
order to produce the assumed rotation-amplitudes. Let 0;=u/, (¢=1, 2,---, m), 
be the rotation-amplitudes and A; the corresponding amplitudes of the couples. 
Furthermore, let B; be the amplitudes of the exterior couples obtained by the first 
step of the method of moment balancing. Then, if the assumed 0; represent the actual 
configuration enforced by the loads B;, A; —B,;=0; but in general 


A; — B; = Ci, (17) 


where C; is the residual moment-amplitude. 
Amongst all possible systems 6; the actual one minimizes the energy function 


7 See, for example, J. I. Parcel and G. A. Maney, Statically indeterminate stresses, John Wiley and 
Sons, 1936, p. 149. 
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@na/PIl’ Oc/PU’ @cs/Pll’ 6cp/PU' 
.000 | 381 — .000 
.307 —.074 
| — .268 
000 | 113 —.074 —.074 
022 | —.091 
| .079 
022 | 022 .005 —.074 
— .064 015 
055 
022 077 — .059 — .059 
O11 —.044 
.038 
033 033 — .021 — .059 
— .020 .008 
.027 
033 060 —.051 —.051 
.005 — .022 
.019 
038 .038 — .032 —.051 
—.015 004 
| 013 
038 | 051 — .047 — .047 
003 — .010 
| .008 
denied | 
041 | 041 — .039 — .047 
| — .006 002 
| .006 
041 047 — .045 — .045 
001 — .005 
.005 
042 | 042 — .040 — .045 
| — .004 001 
| 003 
042 | 045 — .044 — .044 
.001 — .002 
043 — .044 — .044 





64n=Opc=0. 
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Balancing angle changes Balancing moments 


= 
| 
| 








Maz/PI | 116 118 
Mpa/P1 134 135 
Msc/Pi | — .138 — .135 
Mcs/PI .132 135 
Mcp PI | — .136 | — .135 
Moc/Pl —.119 —.118 
H —= ; > 5400. = p BO}. (18) 
i,k=1 k=l 


The first term on the right side represents the internal energy, 


n I l 
4 Zz a;,00;. = {er f (u")*dx — pot f wash (19) 
0 0 


t,k=1 


where the right hand sum is to be taken over all members of the structure. The rela- 
tion (19) arises from the fact that, for any member, u” and u can be expressed linearly 
in terms of the rotation-amplitudes at the ends of this member. Note that aj, =a; 
and A F =). rik Ox. 

Let us denote the first natural frequency of the structure by w:. The values 
, an then can be shown to be positive as long as w <a. Indeed, by Ray- 


9 l l 
wi Ss |x er (u”)? az| /|x+f waz], (20) 
0 0 


where the sums are to be extended over all members of the structure. As the function 
u in (20) let us take the displacements corresponding to 6,;=1, 6.=@;= --- =0,=0. 
From (19) and (20) together with the condition w<w, it is then clear that ay, >0. 
Similarly a2 >0, a33>0, - + +, @an>O0. 

Let a first set of values 6,=6{” be given and compute the corresponding residual 
moment-amplitudes C!”. Suppose the subscripts 1, 2, - - - , 2 to be arranged in such 
a manner that | C{”| =| C{”|, (¢=2, 3, - - - , 2). We now define a second set of values 
6? which differs from the first one only in so far as the value of 6; is concerned: 


O11, Meo, °° ° 


leigh's principle 


: = +¢, i =, , (¢ = 2,3,---,#). 
We propose to determine ¢ in such a manner that the value of H is decreased.* We 
have 
(2) me = (1) : O11 2 () @i1 2 
H(@ )— H(6 ) = | Dandie — Blo+ 6 =Ci ¢+ 7" . (21) 
k=l 
Taking 


* See G. Temple, The General theory of relaxation methods applied to linear systems, Proc. R. Soc. of 
London, Ser. A, 169, 476, (1938-1939). 
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(Q) 
Ci 
o=-- (22) 
E 11 
we obtain 
(2) (1) (1) 42 
He )-—H@ )=—3[Ci | /an (23) 
which is certainly negative as long as w <w;. The residual moment-amplitude C? cor- 


responding to the new values 6!” equals 


(1) 


(2) = (2) 
C. = Dad —Bi=Ci + and = 0. 
k=1 


This shows that the choice of @ according to (22) corresponds precisely to the process 
of moment balancing where in each step the greatest absolute residual moment is 
“liquidated.” For the next step the subscripts i have to be rearranged, so that C® 
is the greatest absolute residual moment. Continuing in this way we obtain a decreas- 
ing sequence of values of 7. If we simplify our notation by writing H® instead of 
I1(@™)), this sequence becomes 


H%> H>--->H>--- > Haw 


with 
(p),2, (rp) 


(p+1) (p) 
H —-H =-—-3(Ci |/an <0. 
Here aj’ has been written instead of the ay of (23), since as a consequence of the re- 
arrangement of the subscripts the value of this quantity changes from step to step. 
Now aj)’ is positive and can assume only a finite number of different values (m at 
the most). Furthermore, the sequence H) is decreasing monotonically and is 
bounded from below by Hmin.. Therefore 


lim [c,”’]? 


po 


= 0. 


Since C;”’ is the greatest absolute residual moment in the pth step, this means that 
ultimately all residual moments will disappear. The structure is then completely bal- 
anced. 

This convergence may be rather slow, especially if w is near w,. For example, com- 
pare the 12 balancings used in Example 2, when A =3.30, to the 5 needed in Example 
1, for the same structure when A =0. For this structure A; =3.55. 

The method of balancing angle changes may not always converge when w <a, 
as will be seen if Example 4 is attempted when A\=3.30. Usually the method of bal- 


ancing moments converges more rapidly than the method of balancing angle changes. 








250 


AN APPLICATION OF THE METHOD OF THE 
ACCELERATION POTENTIAL* 


BY 


J. LEHNER (Cornell University) anp C. MARK (University of Manitoba)** 


In this paper Prandtl’s theory of the acceleration potential is used, in conjunction 
with conformal mapping, in order to determine the pressure distribution on a sym- 
metrical control surface consisting of a fin and flap separated by a gap of finite width, 
under the assumption of a steady irrotational flow of an incompressible perfect fluid. 
The method used is essentially that by which M. A. Biot! recently derived the well- 
known formulae for lift and moment of a symmetrical airfoil with flap in a remarkably 
simple manner. The present problem is considerably more complicated than that of a 
single airfoil; but it is still possible to obtain formulae in closed form. In the case where 
the gap between fin and flap is large in proportion to the chord of either, the formulae 
obtained here for the pressure distribution on the fin, or on the flap, do not differ 
materially from those used for a single symmetrical airfoil. 

For treatments of this or related problems by the classical velocity-potential 
method, see I. Fliigge-Lotz and I. Ginzel, Die ebene Strimung um ein geknicktes Profil 
mit Spailt, Ingenieur-Archiv 11, 268-292 (1940), which also contains references to 
earlier studies. Fliigge-Lotz and Ginzel do not restrict themselves to the symmetrical 
case (fin and flap of equal length) as is done in this paper. The method by which they 
obtain the complex velocity potential is essentially the same as the one we use to 
derive the acceleration potential. They compute the pressure distribution for an un- 
symmetrical split wing rather than the total lift and moment, so that a comparison of 
their numerical results with ours is not practicable. 

A paper by Kutta in the Sitzber. Bayerische Akad. of 1911 considers the special 
case in which the two airfoils have the same angle of attack. Our results agree with 
his if we make the identification sin 2a=2a, where a@ is the angle of attack. 

The first part of the paper gives a description of the methods and results; while 
some details of the mathematical methods used are given in the second part. 


I. GENERAL DESCRIPTION OF METHODS AND RESULTS 


1. The acceleration potential. The equation of motion of an incompressible per- 
fect fluid of density p is 


pad = — grad p, (1.1) 


where @ denotes the acceleration vector, and p the pressure. According to Prandtl,’ 

* Received Jan. 18, 1943. 

** This paper was prepared at the suggestion of Professor W. Prager while the authors were partici- 
pants in the Program of Advanced Instruction and Research in Mechanics at Brown University, Summer 
1942. The authors are greatly indebted to Dr. L. Bers for valuable suggestions. 

1M. A. Biot, Some Simplified Methods in Airfoil Theory, Journal of the Aeronautical Sciences, 9, 
185-190 (1942). 

2 L. Prandtl, Beitrag zur Theorie der tragenden Fldche, Zeitschrift f. angew. Math. u. Mech., 16, 360- 


361 (1936). 
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the function 
¢ = — (1.2) 
p 
may be called the acceleration potential, since d=grad ¢. 

We shall consider the steady irrotational plane flow around a symmetrical control 
surface consisting of a fin and flap separated by a finite gap. Taking the axis of x 
parallel to the velocity U which prevails at an infinite distance from the control 
surface, we write the components of the velocity vector g as 


qz=U+u, dy = 2. (1.3) 


In the case of a thin profile with a small angle of attack, the terms u and v can be as- 
sumed to be small compared with the velocity U of the undisturbed stream. 
For a steady flow, the components of the acceleration vector are 


— Ee ee LE O4y Ody 
@s = 42-—— —-y; = 9g, — seam 
ax ay ee ae a 
Introducing the acceleration potential g on the left sides and the expressions (1.3) on 
the right sides of these relations, and neglecting terms of the second order in u and 2, 
we obtain 
09 Ou 0g Ov 
dane, aid U ia U 


Ox Ox oy Ox 
According to the condition of incompressibility, 0u/dx+0v/dy =0, we have 0u/0x = 
—0uv/dy. Substituting this in the first of our equations gives 





ae av dy v 
—=-U—, —=U—.- (1.4) 
Ox oy oy Ox 
Elimination of v between these two equations leads to 
ao dy 
vam en (1.5) 
ax? ay? 


The acceleration potential ¢ is thus seen to satisfy the Laplace equation and, conse- 
quently, can be taken as the real part of an analytic function, f(x+ty) =¢(x, y) 
+iy(x, y), of the complex variable x+iy. M. A. Biot has shown that the conjugate 
function of the acceleration potential, ¥(x, y), also has an immediate physical 
significance. Indeed, from (1.4) and the well-known Cauchy-Riemann relations, 
d¢/dx =dy/dy, A¢/dy = —AY/dx, it follows that the function y may be defined so 
as to equal — Uv. 

In the rest of this section we shall speak of only one airfoil on the understanding 
that what is said applies equally to the fin and the flap. 

The undisturbed flow U along the x-axis involves a normal velocity + Ua at the 
surface of a thin straight airfoil having the small angle of attack a. The boundary con- 
dition of tangential flow on the surface will be satisfied by a function giving the veloc- 
ity — Ua normal to the surface. As is usual in the linear approximation, this boundary 
condition is applied at the x-axis instead of at the surface; so that we now take a part 
of the x-axis to represent the airfoil (see Fig. 1). We have, then, v= — Ua; and hence, 


y= — Uv= U*a; (1.6) 
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that is, Y is constant on the surface of the airfoil. 

The problem of finding the pressures, and hence the lift on the airfoil, may now 
be restated. One seeks an analytic function of the complex variable x+7y the imagi- 
nary part of which satisfies the condition (1.6) and, since y= — Uv, vanishes at in- 
finity. The real part of this analytic function may then be used as the acceleration 
potential from which the pressures may be obtained by (1.2). Taking into account the 
properties of harmonic functions, the fact that Y must vanish at infinity in the z-plane 
(g=x-+7y) makes it necessary that W have at least one singularity. Aerodynamical 
considerations indicate the leading edge of the airfoil as the obvious location of this 
singularity, partly by the analogy with the classical thin-wing theory in which the 
velocity turned out to be infinite at the leading edge. As in the case studied by Biot, 
we assume this singularity to be a source-sink 
doublet with axis parallel to the y-axis. 

To simplify further the determination of 
y, the segment of the x-axis now used to rep- 
resent the airfoil is mapped conformally into 
ry a circle in the w-plane. 
| 2. The mapping. We treat only the sym- 
metrical case in which the fin and flap are 
of equal length.* Let these be represented by 
slits along the real axis: the fin from —1/k’ 
to —1, the flap from 1 to 1/k’, (0<k’ <1). 
The z-plane exterior to the two slits can be 
mapped into the interior of a circular ring in 
the w-plane, as indicated in Fig. 1. The radius 
of the outer circle may be taken as unity, the 
other radius, R, being then fixed. The slits 
transform into the boundary circles of the 
ring. The function giving the required map- 
ping can be written in closed form using ellip- 


U—. 


-+x 
DCI/k') 


(b): Z-PLANE 








— i 
BCI) 


AGI/k') C1) 


tic functions (see §6). 

In the mapping, the upper edge of the slit 
AB (Fig. 1) goes into the upper semicircle 
A’B’, the lower edge of the slit into the lower 
semicircle. Moreover, points such as G and H 





(c): w-PLANE 


Fic. 1. The actual fin and flap are indi- 


cated in (a); the slits used to represent them 
in (b): it is on these slits that the condition 
(1.6) is actually met. (c) shows the circles into 
which (b) is mapped. It is convenient to use 
temporarily several sets of polar coordinates 
in the w-plane: (r, 0), with origin at center of 
circles; (71, #:), with origin at C’; (r2, 62), with 
origin at A’. 





on AB, G being on the upper edge of the slit 
and H on the lower edge immediately be- 
low G, map into points such as G’(R, @) and 
H'(R, —6) for which the values of @ are equal 
but opposite in sign. Similarly for CD. For 
uniqueness, we define the coordinate @ in the 
w-plane so that | | <r. The point z= 
maps into © '(./R, 7). 


* The unsymmetrical case could be handled just as easily at this stage. Any fin-flap arrangement 





could be mapped into the symmetrical case by a preliminary linear transformation of the type 
z’=(az+b)/(cz+d); subsequent steps being the same as here. It is in meeting the condition ¥(#) =0 
that the symmetrical case is notably easier to treat; unsymmetrical cases might require numerical han- 


dling from that point on. 
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The radius of the inner circle (transform of the fin) is given by 

R = ¢~?*K'/K, (2.1) 

where K and K’ are the complete elliptic integrals of the first kind (modulus &), k be- 

ing the complementary modulus to k’ (k?+k’?=1). The relation between the coordi- 

nate z of a point on one of the slits in Fig. 1(b) and the coordinate 6 of the corresponding 
point in Fig. 1(c) is 

1 1 03,4(30, \/R) 

;  Samnanonees se 2 -) ———e 

dn(K6/z, k) Vk’ 33(30, VR) 


= + 


Z (2.2) 
in which the plus sign holds for z on the flap. The notation of the elliptic and 3-func- 
tions is that used in Whittaker and Watson, Modern Analysis (Cambridge Univ. 
Press, 4th Ed., 1927; Chap. XXI and XXII). 

3. The potential. The problem now is to find a potential which will vanish at in- 
finity and have a singularity of the doublet type at the leading edges of fin and flap 
(the points A’ and C’ in the w-plane), and of which the imaginary part y will satisfy 
condition (1.6) separately for each airfoil; namely, 


y = Ua onr=R, and y = U’8 on sr = 1.* (3.1) 


In addition, corresponding to the Kutta-Joukowski condition in the classical theory, 
the acceleration potential must be continuous at the trailing edges of fin and flap. 

We place a plane doublet of strength a; at C’ and one of strength a2 at A’, the axis 
of each doublet being perpendicular to the real axis. The values of a; and a2 will be 
determined later by (3.1). The complex potential at the point &, in the plane of a com- 
plex variable u, due to a doublet of strength m at the point u=0 whose axis makes 
an angle n with the real axis is 

mé—'e", 


Thus the potential due to the doublet at A’, for example, is given by (see Fig. 1(c)) 
dorz 'e'(*/2—-%) = qorz! sin 02 + taorz! cos Oo. 


The imaginary part of this potential, aarz! cos 02, has the constant value a2/2R on the 
circle r=R, but is not constant on r=1. Similarly, mutatis mutandis, for ayr—' cos 4,. 
The function a,r77! cos 6:+aerz! cos 62 will not, then, serve as the ¥ function, but it 
may be modified as follows. 

We obtain the harmonic function F(r, 6) which, to within an additive constant, 
takes the values —a,r7! cos 0;0n r=R, and —aerz' cos 62 on r=1. Then the function 


Y = a,rz' cos 0; + aorz' cos 02 + F(r, 0) + ¢, (3.2) 
where aj, a, and the available constant ¢ are chosen so that (3.1) and 
¥(2’) = W(WV/R, x) = 0 (3.3) 


are satisfied, will be the imaginary part of the potential required. 

Equations (3.1) and (3.3) are satisfied (see §7) when 

* When (r, @) is the map of z we shall use ¥(z) and y(r, @) indiscriminately where no misunderstanding 
can arise. The variable indicates the plane (and the point in the plane) at which the y function is to be 
considered. Similarly for other functions with a physical significance. 
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kK 3 ade K 1 
u(t )+3(E-2) 
2a 4 R \2r 4 


fo 








2K(a — £) — 
a, = y(=S—* mete) (3.4) 
K(1 + #’) 
a, a — 8) + x(a + s 
R K(1 + ’) 


(The apparent difference between the ways in which a; and az appear in (3.4) is due 
to the fact that a2 is on a circle of radius R.) 
The function F(r, @) in (3.2) was obtained in the form of a Fourier series: 


F(r, 0) = : (A,r”™ + B,r-") cos né. 
1 


The conjugate harmonic function may then be written immediately as: 


— > (Ayr™ — Byr-*) sin nd. 


1 
The conjugate function to 
ayrz' cos 6; + aerz! cos 62 is — ayrz) sin 0; + aorz! sin Bo, 
(the minus sign being due to the fact that 6,1, as defined in Fig. 1(c), is measured in the 
negative sense). In this way the acceleration potential is obtained to within an addi- 
tive constant b (which disappears in evaluating the lift). Dropping the auxiliary co- 


ordinates, we have finally (see §7) 





he fi ae pepe. os tani % os R) +5 
on the n: p\ = —— (89, ns. a -omimmne (3 : : 
? 2R v2 . 2 v4 P 
(3.5) 
the fla (1, 6) ah ap 2 Of us RV +5 
on e : D , SS =e, — — — (490, ; 
J 4 2 vy , 2R v3 , 


in which #/ /3;(36, R) is written for 


1 d 
ae Swe 
3:(40,.R) d(}0) ? 





4. The lift and moment. From (1.2) we have p= —pg. The lift on the air-foil at 
such a point as G is given by 


I(z) = p(H) — pG) = p(B’) — pC) = ple) — o(— ®)] = 2p(y(6) — ) = 1). 


Thus, on the fin, 


I(R, 8) ( a: 84 a, R) ve (40 R)) (4.1) 
’ = p 7 vere . 2" 9 sete, | eens ’ ’ ° 
Rw” ‘— 
on the flap, 
oy ; a2 og 
(1, 0) = p( — a, — (de, 8) - = = (o, 2); (4.2) 
vy R #3 
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in which a; and az are given by (3.4). 
The total lift (on the flap, say) would then be given by 


1/k’ r dz 
a -{ \(2)de = f 1(0) — do. (4.3) 
1 0 dé 


The moment of the flap about the point z=0 is 
1/k? 
M (0) -{ l(z)zdz; (4.4) 
1 


and about any point 2, it is 


1/k’ 
M2(z1) -f l(z)(z2 — 2,;)dz = M2(0) — 2,L>. (4.5) 
1 


The functions involved in these integrals are found in (4.2) and (2.2). It will be ob- 
served that the J-functions appearing have different parameters. 

5. Results. The above integrations can be carried out exactly. The method is de- 
scribed in §8, and the results are: 


a 6 
p1U~*L2=— [x(1— hk’) +2Kk’ — E] + —— (1—’)(4?+4K?'), 
k’ 2rk’ 7 
a . 
o1U-2M (0) = ro (1—k’)[2E— $2 (1—#’)] 


+ [r(1—k’) { Kk’ +E—i x (1—k’)} +(Kk’—E)?*], (5.1) 


rk’? 


p'U-*L, = - [2(E— Kk’) —(1—#’)(2Kk’—7) ] 


6 
ee [w{2E—K(1+k’2)} —4(1—’)(2Kk’—x)(2K+2)], 
where L2, M2(0) have the meanings given in §4, Z; is the total lift on the fin, E is the 
complete elliptic integral of the second kind modulus k, and 6=8—a. 

As a limiting case we consider the situation where the gap between fin and flap 
is so large in proportion to the chord of the fin (or flap) that they may be expected to 
act as independent airfoils. In this case k’=1, k=0. Reference to tables of elliptic 
integrals (e.g., Jahnke, Emde, Tables of Functions, Teubner, 2nd Ed., 1938) or direct 
integration gives K =E=7/2;so that, from (2.1) we have R=0. As an approximation, 
we neglect powers of R, and our formulae become (Whittaker, Watson, p. 489) 


¢,= —28U?, a;/R = 2aU’, 


1(1, 6) = 2BpU? cot 30. (5.2) 


The expression in (5.2) is the usual one for the lift on a single straight wing with 
angle of attack 8. For example: when the gap between fin and flap equals 8 times the 
chord of either, R=0.0008; and evaluating (4.2) we get, in this case, 











I(1, 0) = pU2(2.1 8 — 0.1 a) cot 48, 
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(5.3) 


which differs only slightly from (5.2) unless |B—a| is large; while for B=a, (5.3) re- 


duces to (5.2): 


Some numerical examples are given in the table which follows. Cz, Cu, Cp are the 





ke 


/ 





Pp 





usual lift and moment coefficients: lift coefficient, 
C,=L/p(U?/2)c, ¢ being the total chord from 
leading edge of fin to trailing edge of flap; mo- 
ment coefficient of the flap about its leading edge, 
Cu = M2(1)/p(U?/2)c?; moment coefficient of flap 
about P, Cp = Cy—Cz,HP/c where HP is the dis- 














Fic. 2. 


tance from the leading edge of the flap to P (see 
Fig. 2). The subscripts 1 and 2 indicate fin and 
flap, respectively, while L indicates lift on either. 


In the table HP is expressed as a fraction of c’, the chord of the flap; and 6=B—a 


and a@ are expressed in radians. 





TABLE 


(a and 6 in radians) 








0.3 




















s/e’ 0.1 0.2 0.4 0.5 
: Lirt COEFFICIENTS 

Fin: Cz, | 4.570 +2.505 | 4.19a@ +2.055 | 3.88a +1.676 3.63a +1.426 3.4la +1.236 
Flap: Cz, 11.41a~@ (2.056 | 1.51a 4+2.185 | 1.58a +2.226 1.6la +2.215 |1.62a +2.185 

HP Moment COEFFICIENTS 
—— | — Se a = z Se caso — = — 
0; Cy | 0.210 +0.285 | 0.20a@ +0.265 | 0.192 +0.255 | 0.18@ 40.245 |0.17a +0.225 
O.1¢': C.s 0.150 +0.185 | 0.13@ +0.166 | 0.12a@ +0.156 0.12~@ +0.145 (|0.1la +0.146 
0.26’: Cs 0.080a+0.08065 | 0.065a+0.0646 | 0.056a+0.0585 | 0.053a+0.0525 0.044a-+0.0495 
@.25e": Cs 0.046a+0.03265 | 0.031a+0.0156 | 0.022a+0.00945 | 0.017a+0.00606 | 0.012a+0.004865 

i = | 





Finally, it might be of interest to note that one can locate a point about which the 
moment on the flap would be proportional to 6. Although, if hinged at this point, the 
flap would not remain in the position 6=0 without some restraint, still, when in this 


position, the moment on it would be 


s/c’ = 0.5, 
s/c’ = 0.1, 
II. 


zero. Thus if: 
— 0.0146; 
— 0.0366. 


0.270c’, 


0.313c’, 


C .270 


C313 


HP 
HP = 


MATHEMATICAL APPENDIX 


The following references will be used: 
W =Whittaker and Watson, Modern Analysis, Cambridge, University Press, 4 ed., 


1927. 


J =Jahnke and Emde, Tables of Functions, Leipzig, Teubner, 3 ed., 1938. 
6. The mapping. (See §2.) The mapping is accomplished in two steps. First, we 


map the slit z-plane into the interior of a rectangle in the {-plane by the function 
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= | i [(1 — 2*)(1 — #2") }-V8de. (6.1) 
0 


This rectangle has vertices at +K’(k) +7K(&), the elliptic integrals being expressed as 
functions of the modulus k complementary to k’. The correspondence is: z=1—>¢ = K’, 
z=k’-!—¢ =K’+iK, (plus, for the upper boundary of the cut), etc. The integral be- 
ing multiply valued, the mapping is repeated infinitely often’ in the {-plane, covering 
the plane without gaps by a net of congruent rectangles. 

The second step of the mapping is 


w = exp [x(¢ — K’)/K]. (6. 2) 


This is periodic with period 27K, hence w takes both horizontal sides of the rectangle 
into the part of the real w-axis between —1 and —R, (R=exp (—27K’/K)), i.e., into 
the segment D’A’ (Fig. 1(c)). The vertical sides are mapped into the circles | | =R, 
| w| =1. Because of the periodicity of the exponential function, the full network of 
rectangles in the ¢-plane is mapped into the ring RS|w]| S1, congruent points of 
different rectangles going into the same point of the ring. Thus the mapping between 
the z- and w-planes is one-to-one. 
The inverse mapping is, with w=re”, (W, p. 492), 


K 
z = sn ({, k’) = sn (~ log w + K’, ¥’) = sn (— (10 + log r) + K’, #). 
T T 


For the two important cases, r=1, r=R, the last formula reduces to (2.2) if we trans- 
form the elliptic function to one having a real argument (W, pp. 500-506). Note that 
when the elliptic function is expressed in terms of 3-functions, the parameter, q, of 
the 3-functions is exp (—7K’/K) = R!/? (W, p. 479, ex. 3). 

7. The potential. (See §3.) The function F(r, 0), assumed as a Fourier series 


F(r, 0) = >> (Ayr™ + Bar-") cos n8, (7.1) 
1 


must, to within an additive constant, satisfy the boundary conditions 
> a,(1 — R cos 6) 
1+ R?—2Rcos@ 
— a2(R + cos 6) 
1+ R?+ 2R cos 





F(R, 0) = — ayrz! cos 6; = 
(7.2) 





F(1, 0) = — aerz! cos 02 = 


the expressions on the right being obtained from Fig. 1(c). But 





1 — Rcos #6 1 = 
- — = R ———— = } R* cos nd, 
1+ R? — 2R cos 0 1 — Re*® oat 
R + cos 0 e~# * 
——— = { ——— — = )) (— 1)"R" cos (n + 1)8, 
1+ R?+ 2R cos 0 1+ Re-*® m4 


where § denotes the real part; and hence, substituting the last two expressions in 
(7.2) and comparing with (7.1), we get 


A,R" + B,R-*® = — a,R", A, + B, = — ao(— 1)" R™, (n 2 1), (7.3) 
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with solutions 


a,R*" + a(— 1)" 


: a,R?" —+ ae — 1)*R**-! 
1 — R™ ° 


- —- — —- (7.4) 
1— R2" 


e229 





On substituting these values for A,, Bn, or, more conveniently Eqs. (7.3), in the 
definition (7.1) of F(r, 6), we find 


F(R, 0) = — ary! cos 0; + ay, F(1, 0) = — derz! cos 02; (7.5) 


so that F(r, @) satisfies the required conditions. The constant c in (3.2) must be chosen 
so that y vanishes at the infinite point of the z-plane; i.e., (see (3.3)), 


y(R! _ ) = [— ari 1 cos 6; + dor! cos ee 
+ © (— 1)"(A,RY2" + BR") + ¢ = 0. 
1 


The value of the bracket is obtained from Fig. 1(c); and using (7.4) we find 


, oo R32 ts R") 
fe ee een eee te (= 1) eee 
1+ RU2 " RV — RV) md 1 — R* 


a, 2 8/1 — R*) 


Rj 1 — R* 





(7.6) 





For the evaluation of the infinite sums in (7.6) we need certain formulae from the 
theory of elliptic functions. These are 
=) q” K 1 ri 
a ————$—— = — — — (W, p. 511, ex. 1), 
rT itg¢g* 2e 4 
ca — { nan K R’ 1 
> i J I ad (ibid., ex. 2). 
1 1 + ” aed 2r 4 


Using R=q? (cf. end of §6), and making some obvious algebraic reductions, we find 
for c the value given in (3.4). 

The values of a; and a2 are found from (3.1) which with (3.2) and (7.5) give the 
equations 


$d2 R+ate = U?a, taun+c= UB. 


The solutions of these are contained in (3.4); and thus y is completely determined. 
The acceleration potential g is now obtained as the harmonic conjugate to the 
function ¥. We have, (see §3) 


oo 
g = — aj sin 0; + aerz sin 02 — >, (Anr® — Byr—") sin no + 6, 
1 
where 3 is an arbitrary constant; or, in terms of r and 8, 
— ayr sin 0 der sin @ 
olf, 0) =- nt EE Soe 


1+72—2rcos6 | R?+ 42+ 2rR cos 8 . 
(77) 


— > (A,r" — B,r-*) sin n0 + b. 
1 
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For the purpose of evaluating the lift we need only g(R, 6) and ¢(1, 6). These are given 
in (3.5); we establish the expression for g(1, 0), that for g(R, 6) being obtained in a 
similar way. 

If, in (7.7), we set r=1, use (7.4), and collect terms we find 








ay oo R2" 
o(1, 0) = — — <cot 30 + 4 ——— sin no 
, 2 \ . ms 1 — R 
a2 — Rsin#@ * : 
_-— \ - _ pS (— 1)"R* sin n6 
R (i + R? + 2R cos 6 1 
2 (— 1)"R" 
+2>)> — =. sin no + b. (7.8) 
1 1 — R** 
But 
— Rsin#@ Re-*® as 
—————— = $ (— ) = Z (— 1)"R* sin n8, 
1+ R?+ 2R cos 0 1+ Re-*® 1 


where % denotes the imaginary part. Using this in (7.8) we obtain 
a a) R2" 
g(1,0) = — — {cot 49+ 4 ———— sin no 
2 ' ~ 1— R* 


ao = (— 1)"R” \ 
—— <4 ———— sin n> + b. 
atte 


(7.9) 


Finally, the terms in the brackets in (7.9) may be expressed by means of 3-functions. 
Making use of the results of W, p. 489, ex. 12, we get 
ay vi a2 v3 
(1, 0) = —— —— (30, R) — —— — (30, R) + 4, (7.10) 
, ie Re 
as stated in (3.5). It is worth pointing out that the parameter of the d-functions met 
here is R, whereas in equation (2.2) the parameter is R*’?. 
8. Some definite integrals. (See §4.) The integrals in §4 aie combinations of the 
following eight integrals: 


¥ o,, dz 
a -{ “" (19, R) — dé, (8.11) 
0 Gus dé 
a dz 
Tus -f — (36, R)z — dd, (ms = 1, 2, 3, 4), (8.12) 
0 m dé 


where z and @ are connected by (2.2). For their evaluation we can make use of the 


following formulae: 


1 od Az 

a aoe (40, R) = ——— + B, (8.21) 
d@ vy z—1 

d v7 A 

— — (i, R) = ——— + B, (8.22) 


dé ve 1 — k'z 
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dee R) : +B 
— — (20, = ——__—_ ; 8.23 
do 3, ° 1+ k's ieee 

d v4 (40 R) Az 4 I 
— — (10, R) = ——— B, 8.24 
de 4 } Z + 1 \ ) 

where 
1 
z= — (8.25) 
dn( (K6, Tr, r, k) 

A = — K2k?/r?, B = K(K — E)/z’. (8.26) 


We shall prove only (8.21), the others being established in exactly the same way. 

From the quasi-periodic properties of the J-functions (W, p. 465, ex. 4) we see 
that the left member of (8.21) is a doubly periodic function with periods 27, 277, 
where 7 =27K’/K. We may restrict ourselves to a single period rectangle, say the one 
having the origin as the southwest vertex. Here the left member of (8.21) is regular 
except for a pole with principal part —2/6? (W, p. 466, p. 489, ex. 12, and the Laurent 
series for the cotangent). But, using (8.25) and W, p. 504. 


k? K6? es 
pe er) ae Re: 


2 x’? 


(8.3) 
Az 2An? 4" 2 
anon tn semereten  AQOY) ee — + Ctmat, + - >>. 
s—1 kK? @° 
Furthermore, Az/(z—1) is regular at points in the period rectangle other than @=0, 
for z is an elliptic function of order 2 and (8.3) shows that it takes the value 1 twice 
at @=0. 
We see, therefore, that the difference 


is a doubly periodic function without singularities; thus it is a constant, B, which we 
evaluate at @=7. Using W, p. 489, ex. 12, we have 


d vi i 2 ) , y+ar’ nR* ” cos nO 
— —— (30, KR) = — § CSc’ § peer 
dé vy we r 1 1 — R?* 
and at 6=7, 
d 2 (-— 1) “nR* 
— — (je, 2) = —34+4>, ———_ 
a X 1 — Re 
7 - 1)"ng*" =. (— 1)*sq** 
ap 42f DF COW pn 
1 i—q* 1 1 + g** 


The sums are evaluated by reference to W, p. 535, ex. 57, and p. 512, second formula 
differentiated. We obtain 


d df K*k’ KE 
— — (4x, R) = —-——- —.- (8.4) 
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Also, by (8.25), z=k’—! where 0=7. Hence, 


K K 
B= — [-— Ki — E+ K#k'-(k’ — 11] = — (K - B), (k? + k”? = 1), 
= 


Tv 
as in (8.26). 
As before, we shall evaluate only one of the integrals, say J; in (8.11). This is ac- 
complished by an integration by parts. 


oy dz 
I; -f —_ (16, R) — de 
0 Vv dé 


Il 


oy . *d oy 
(10, R)z(6)| — f — — (30, R)z(6)d0 
vy lo 0 dé vy; 


oy uk od wy dé 
— | — (10, R) -f — — (10, R)z — ds. (8.5) 
vy 60 1 dé vy dz 


In the extreme right member of (8.5) both the integrated part and the integral are 
infinite, but their difference, considered as the limit 


oy sky od Bf dé 
im ( - (30, R) -{ — — (30, R)z — dz}, 
«0 vy Ome z(e) dé vy dz 

is finite. 


The integration can now be carried out without difficulty. In the integral in the 
right member of (8.5) use is made of (8.21), while d0/dz is obtained from (8.25), the 
result being then expressed as a function of z. It is necessary to show that the singular 
contributions from the two terms of the right member of (8.5) cancel. This is done 
conveniently by employing their Laurent expansions. The integrals which present 
themselves are at worst elliptic integrals of the first two kinds, and can be found in J, 
pp. 52-56. In this way we obtain the results of §5. 
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PERIODIC PROPERTIES OF THE SEMI-PERMANENT 
ATMOSPHERIC PRESSURE SYSTEMS* 


BY 
H. J. STEWART 
California Institute of Technology 


The outstanding features of the general circulation of the atmosphere are the 
belts of westerly winds and, on the equatorial side of these, the system of semi- 
permanent sub-tropical high pressure areas. In a previous paper! the author has dis- 
cussed the problem of the formation of such high pressure systems. In particular, it 
was shown that these systems probably represent dynamically stable concentrations 
of vorticity similar to the Karman “vortex street” which is formed behind any two- 
dimensional bluff body over a wide range of values of the Reynolds number. It now 
appears that a further examination of the periods of the characteristic oscillations of 
such systems is of considerable interest. It is seen that the period of these oscillations 
is of the order of magnitude of years. This indicates that oscillations of this type may 
be of importance in the calculation of the long period displacements of the Pacific or 
Azores high pressure systems. 

It is believed that this is the first time that atmospheric motions have been dis- 
cussed which have a period of the order of magnitude of, but different from, a year. 
Since the weather shows large variations from one year to another, it is apparent that 
such motions must exist; and, since the non-seasonal variation of the only external 
parameter, the solar energy input, is very small, these long period motions must be 
explainable in terms of the free oscillations of the earth’s atmosphere. 

It seems that the horizontal field of motion is of primary importance in determin- 
ing the motion of these large scale systems; so it is assumed that the atmosphere can 
be treated as a single layer of fluid of constant density with the vertical velocities 
being of small importance so that the pressure can be determined from the hydro- 
static equation. It is also assumed that the apparent acceleration is negligible when 
compared to the Coriolis acceleration. In addition the effects of friction and of the 
variation of the Coriolis parameter with latitude are neglected. This latter factor 
means that the fluid motions considered are those taking place on a rotating disc 
rather than on a rotating sphere. 

The notation used in the discussion is as follows: 

x, y = Cartesian coordinates on a rotating disc, u = velocity in x direction, v = veloc- 
ity in y direction, w=angular velocity of the disc, h=depth of the fluid, g=accelera- 
tion due to gravity. 

If the motion could have been started from rest with a uniform depth ho, the 
principle of conservation of the absolute vorticity states that 

Ov Ou 2w 
—— — = — (kh — hy). (1) 
Ox dy ho 


* Received August 4, 1943. 
1 Stewart, H. J., Proc. Nat. Acad. of Sci., 26, 604 (1940). 
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If the velocity components are eliminated from this equation by means of the geo- 
strophic wind equations, 


Oh Oh 
— 2w= —g—) 2wu= — g—, (2) 
Ox oy 


an expression determining the depth of the atmosphere (i.e. the sea-level pressure) 
is obtained. This is 

Oh Oh 4? 
— — —(h— hy) = 0. (3) 


Ox? dy? = gho 


This equation can be further simplified by the introduction of dimensionless variables, 
X =2wx/4/ghy, Y=2wy/s/gh) and n=(h—ho)/ho. With these new variables, Eq. (3) 


becomes 
0°n 0°” 


~gme 4 
ox? ey? © 





In terms of the dimensionless depth and dimensionless velocities defined by 
U=u// gh, and V=0/1/gho, the geostrophic wind equation can be rewritten as 
On On 
—=YF, —=-U. (5) 
ox oY 
The only solution of Eq. (4) which vanishes at infinity and which represents flow 
in circles about the origin and thus corresponds to a simple vortex is 


n = aK9(r) (6) 


where a is an arbitrary constant, r=+/X?+ Y? and K,(r) is a modified Bessel function? 
of the second kind. If @ is positive the motion is anticyclonic; if a@ is negative the 
motion is cyclonic. All of the motions considered in the present investigation are built 
up through superposition of vortices of this type. From Eq. (5), this vortex has a di- 
mensionless tangential velocity, us, given by 

on 

u = — aK,(r) = —- (7) 

or 
The geostrophic wind equations used in the above development can be shown to be 
valid unless r<1. 

Based on a homogeneous atmosphere having a mean sea-level pressure and density 
of 1.013 10° dynes/cm? and 1.22X10-* gm/cm* respectively, the same as the stand- 
ard atmosphere, the characteristic velocities and distances used above to produce di- 
mensionless variables are \/ gh, = 2.87 X 104 cm/sec and +/ gh,/2w =1.97 X108 cm. Ata 
distance of 2000 km. from the center of the Pacific or Azores high pressure systems, 
the characteristic velocity is of the order of 10 meters/second. Since K,(1) =0.602, 
this indicates that these anticyclones have a strength such that @ is approximately 
0.06. 

If the interaction between the northern and southern hemispheres is neglected, 
the ring of subtropical anticyclones can be roughly represented by N equal anti- 
cyclones of the type given by Eq. (5) which are placed on a ring of radius a and spaced 





2 Grey, Mathews and MacRobert, Bessel Functions, Macmillan and Co., London, 1931. 
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at equal angles t=22/N as shown in Fig. 1. In the northern hemisphere, the Pacific 
and Azores high pressure regions are well defined. They are about 120° of longitude 
apart. There is some evidence of a third such system over India and equidistant from 
the other two; however this evidence is far from conclusive due to the low level inter- 
ference from the monsoon. In the southern hemisphere there are also three such sys- 
tems. The best model is thus obtained with N=3. 
The surface deflection for such a system in its equilibrium state is 
N 


n = a), Kola? + r? — 2ar cos (6 — nr) |*/. (8) 
n=1 
From Eq. 5 the dimensionless velocities in the radial and . 
tangential directions, u, and us respectively, are given 


b 
' : «@ 0 ” 
on on \ 
“h= —-— —» “= —: (9) 
r 06 Or O 
The velocity of any vortex is the velocity at that point a 


due to all of the remaining vortices. From the second of < 
the expressions of Eq. 9, the system shown in Fig. 1 is cE 


seen to have a dimensionless angular velocity Q2 given 
Fic. 1, Ring system of 


by x 
anticyciones. 
a N=! 
Q = —— >> K,[2a sin 4nr] sin 4nr. (10) 
a n=1 


In Table 1 are given the values of 2 and of aQy/ghy, the linear velocity of the vortices, 
for a=0.06, N=3, and a=3.0 and 
; : 3.5. The values of a=3.0 and 3.5 
ANGULAR VELOCITIES OF VORTEX SYSTEMS, : . . 
006. Nut correspond to the ring of subtropi- 

a=v, , «VV Se . . 
cal anticyclones being placed at 


TABLE 1, 





‘ | 2 | Ova latitude 37° and 29° respectively. 
EE ee —_—— — From this result it is seen that such 
3.0 — .000112 | —9.6 cm/sec a vortex system would have a slow 
3.5 — -000037 | —3.7 cm/sec precession to the west. In the at- 


enna oyidal ; 7 mosphere, there is also a region of 
distributed cyclonic vorticity to the north of the westerly winds. It is easily seen that 
this cyclonic vorticity tends to produce an eastward displacement of the subtropical 
highs. It appears that these two displacements cancel one another so that the sys- 
tems are practically stationary, a condition which is undoubtedly also imposed by 
the thermodynamic and topographical factors. Since the westward drift shown in 
Table 1 is very small, no attempt will be made to correct the vortex system shown 
in Fig. 1 in order to take into account the polar cyclonic vorticity. 

If the mth vortex is displaced by a distance Ar, in the radial direction and by an 
angle A@, in the tangential direction, the surface deflection in the displaced condi- 
tion is a 

n = a>, Ko[(a + Arn)? + r? — 2(a + Ar,)r cos (8 — mr — AO,) |*/?. (11) 


n=1 


1943] PERIODIC PROPERTIES OF ATMOSPHERIC PRESSURE SYSTEMS 265 


The velocities of the vortices may be calculated as before by Ea. (9). If the displace- 
ments are small, the changes in velocity of the Nth vortex from the equilibrium value 


indicated by Eq. (10) may be written as 


d N-1 
Au, = r (Ary) = 3 ~ & AraKo(Rn) sin (nr) 
dat n=1 


N-1 K (R, 
+ a> a(Ad,, — Ady) ‘ = 5 ) + Ko(R,) cos? Jnr 
n n (12, 





d N= Ki(R,) 
Aup = a (A@v) + QAry = aAry >> 4 Ko(R,.) sin? 4nr — — cos (nr) 
at n=1 n 


N-1 R, 
+f- a>, ar \K o(R,,) sin? ur + - —_ - ‘\ + 3aA9,Ko(R,) sin om), 


l n 


where 


R,, = 2a sin 3nr. 
[It should be noted that the ¢ in Eq. (12) is a dimensionless time. If the actual time is 


t*, then 

t= 2wf*. (13) 
Expressions similar to Eq. (12) for the velocities of the other vortices could be written 
from symmetry. These would form a set of simultaneous differential equations for the 


displacements. 
If the N equations in each of the two sets indicated in Eq. (12) are added, it is seen 


that 
nad A n -_ 
dt = {2 : } =o 
dit NX , N N-1 1 
) ys a0, = = a: i Lan! ‘ ye | 2x 0(R,) sin? nr + — K,(R,) sin snr 
1 


at l n=1 a 


(14) 


From this it is seen that if a mean value of a is chosen so that } aS Ar, = 0 initially, 
then from Eq. (14) both }°*_,Ar, and >-*_,A6, will remain constant. These results 
correspond to similar equations for two dimensional line vortices which state that 
the impulse of a system having no external forces remains constant.? © 

The disturbed motion of the vortex system as described in Eq. (12) can best be 
discussed by considering the normal modes of oscillation. From the symmetry condi- 
tions, the displacements in each normal mode must be of the form Ar, =Arye™?* and 
Aé,, =AOyve'"* where ¢ characterizes the normal mode and is a member of the series 


2r/N, 4r/N,-+-+, 2n(N—1)/N, 2x. With this notation, Eq. (12) can be written as 
1 d ad “i 
- — (Ary) = AAry — Baddy, §— —(AOw) = CAry + Aadby, (15) 
a dt a dt 

where 


* Lamb, H., Hydrodynamics, 6th edition, Cambridge University Press, London, 1932, 220. 
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1 N—1 
A = 7 :e Ko(R,,)e'"* sin (nr) 


n=1 
W=1 K (R,) 
B= > (1 —ei**) f ae + Ko(R,) cos? oe | (16) 
a UR, f 
+ = > \ , oa K,(R,) * 
C= 304 K0(R,)(1 + ei”*) sin? Ant + a (1 — 2 cos (mr) + e*™*) >. 
n=1 n 


If the amplitudes Ary and A@y are assumed to vary like e‘?‘ and p is the dimen- 
siortless normal frequency, then from Eq. (15), 


p=al— id + V/BC}. (17) 


From Eq. (16), it can be seen that for g=27, A =B=0 and there are thus two zero 
normal frequencies. These two zero frequencies are those shown in Eq. (14). It can 
also be seen that for the specified values of g, A is always a purely imaginary quantity, 
B is always real and not less than zero. From Eq. (17), the condition that the frequen- 
cies be real is that C be real and non-negative. Complex frequencies, of course, charac- 
terize systems in which the amplitudes increase with time and are thus unstable. Now 
C is always real and is always positive for VN <7. If N=7, Cis always positive if a>71. 
For N>7, C is negative for one or more of the given values of g. A value of a>71 
corresponds either to disturbances of such great wave lengths or to motions of such a 
shallow layer of air in the earth’s atmosphere that it probably is of no significance. 
The vortex system is thus stable if N is less than or equal to six. 

The frequencies and modes of oscillation will be discussed in some detail for the 
cases where NV =3 and a=3.0 and 3.5. From Eq. (13) it may be seen that the period 
of an oscillation is given by 


1 
T = = | sidereal days. (18) 


The normal mode of oscillation, from Eq. (15), is given by 

. aé,, addy A — ip/a 

i~ wee. (19) 
Arn Arn B 


The results for a=3.0 and 3.5 are given in Tables 2 and 3, respectively. The two nor- 
mal modes thus show a short period (2000 days to 5000 days) and a long period oscil- 
lation (70,000 days to 250,000 days). The path of the vortex is in each case an ellipse. 
For the short period oscillation the vortices are east of their mean position when 
traveling south and west when traveling north. For the long period oscillation the 
sense of the rotation is reversed. 

Conclusion. The present calculations cannot be considered as a quantitative theory 
of the oscillations of the semi-permanent high pressure systems; they must be con- 
sidered rather as an existence proof. Since the essential features of the model, vorticity 
concentrations at distances of roughly 10,000 km., are also found in the atmosphere, 
motions of this type must exist in the atmosphere. It might be expected that the 
effects of coupling between the systems of the Northern and Southern Hemispheres 
and of any cyclonic vorticity concentrations on the polar sides of the westerly winds 
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TABLE 2, 


NORMAL MobDEs OF OSCILLATION FOR N=3, a=3.0, a=0.06. 
























































2x/3 4n/3 
: 0.002248 —0 .002248 
oom —ti(‘<i«é‘«‘wS*~*~*~*«« 
0.001182 a en 9.001182 
0.000268 | 0.000002 nt 9.000268 “ —0.000002 
1865 | 2.5 xt08 1865 | 2.5x108 
—0.53i | +0 .53% +0 .53i —0.53i 
| 
TABLE 3. 
NorMAL Monks OF OSCILLATION FOR N=3, a=3.5, a=0.06. 
2x/3 4n/3 i pis 
0.00859 =—Sti(<té‘q;~~ (0.000859 
0.001467 el amet ti 0.001467 
0.000375 | 9.900375 <a. 
0.000096 0.000007 0.000096 9.000007 
5,200 7.1108 - . 5,200 74 7.1106 
—0.51i 0.51% 0.51% —0.51% 








would be to decrease the period of the shortest oscillation and to introduce additional 
natural frequencies. No attempt has as yet been made to estimate the magnitude of 
these effects or of the errors involved in using velocity distributions corresponding to 
vortices on a rotating disc rather than to vortices on a rotating sphere and in neglect- 
ing the seasonal variations in the strength of the semi-permanent high pressure sys- 
tems. It is suggested that the present calculations may prove useful as a guide for 
the statistical analysis of empirically obtained data. 
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—NOTES— 


ON HERZBERGER’S DIRECT METHOD IN 
GEOMETRICAL OPTICS* 


By J. L. SYNGE (Ohio State University) 


1. Introduction. In recent papers M. Herzberger'? has developed a “direct 
method” for analytical ray-tracing through an instrument of revolution. At the end 
of the first paper he refers to Hamilton’s method, which he says “leads to an elimina- 
tion problem, hitherto unsolved.” Nevertheless the question arises: What is the 
connection between Herzberger’s approach and that of Hamilton? This question 
is best answered by attacking Herzberger’s problem by the method of Hamilton. 
As we shall see, this is quite feasible. Indeed, if we combine Herzberger’s “direct 
method” with Hamilton’s character function we obtain a very powerful technique. 

Section 2 contains the formulation of the problem of determining the Herzberger 
transformation when Hamilton’s angle-characteristic is known for the instrument in 
question. Herzberger’s identity (AD—BC=1) is obtained immediately. 

In Section 3 the case of a single surface (refracting or reflecting) is considered. It 
is found that the coefficients are connected by a new relation. 

In Section 4 I show how the problem of the sphere may be treated, Herzberger’s 
geometrical approach being replaced by a more systematic analytical method. 

2. The Herzberger transformation. To facilitate comparison with Herzberger's 
work, I shall use his notation. The following table shows the correspondence between 


the notations of Herzberger and Hamilton: 


Herzberger Hamilton 
Coordinates of point on incident rays - a 7,8 
Components of incident ray én ¢ o’, 7’, vu’ 
Coordinates of point on final ray i ef x,y, & 
Component of final ray er ee , %.1 


According to the method of Hamilton there exists an angle-characteristic 7, a 





function of £, 7, &’, n’, such that the equations of the incident and final rays are® 
- aa , Tel [sel Gh 
r— s/f = Ty, g—st/fi =—Ty, 
* ee << (2.1) 
— 2n/f = T,, n/t =—Ty. 


The subscripts denote partial derivatives. 
instrument is of revolution and that the axes Oz, O’z’ lie 


Now suppose that the 


along its axis. Then 7 is a function of the quantities 


* Received May 21, 1943 
1M. Herzberger, Trans. 


Amer. Math. Soc. 53, 


218-229 (1943). 


2M. Herzberger, Quarterly of Applied Mathematics, 1, 69-77 (1943). 
3 J. L. Synge, Geometrical optics, Cambridge, 1937, p. 31. 
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us = (+7), w=’ +m’, us = 3(¢% +77). (2.2) 

Let us write 07 /0u;=T3, etc. Then, by (2.1), the intersections of the rays with the 
planes z=0, 2’ =0, satisfy 

Tst + Tek’, s=u—Té— Tes, 

Tn t+ Ta’, = y = — Tan — Tor’. 


Il 


x 


(2.3) 


Il 


sf 
These equations involve the eight quantities 
, , ¢/ / 


eR: a ae he £59, 2%. 


The basis of the Herzberger method is to express the first set in terms of the second 
set. To do this, we introduce 


u, = 3(2? + y?), Ue = x + yn. (2.4) 
Let us multiply the x, y equations in (2.3) by &, 7, respectively, and add; this gives 
us = 2T3u3 + T4utg. (2.5) 

Rearranging the x, y equations in (2.3), squaring and adding, we get 
Tins = 0, — Tote + Tome. (2.6) 


Supposing JT known as a function of u3, us, us, we have in (2.5), (2.6) two equations 
to determine u,, us in terms of 1%, U2, U3; Suppose the solutions are 


Us = f(u1, U2, Us), Us = g(t, U2, Us). (2.7) 


Making this substitution, we may express 73, 74, 7; as functions of 1, ue, Us. 
Now let us rearrange (2.3) into the Herzberger form: 


x’ = Ax + Bé, t’ = Cx+ Di, (2.8) 
y =Ayt+ Bn, a =Cy+Dn. 
The coefficients are as follows: 
A = — 17,77", B = T,1T,77 — T,, C=Tr, D=-—T;Tz. (2.9) 
We immediately deduce Herzberger’s identity 
AD — BC = 1. (2.10) 


To sum up: Given the angle-characteristic T(us, Us, Us) of an instrument of revolution, 
we obtain the coefficients A, B, C, D of the Herzberger transformation in two steps: 
(i) We solve (2.5), (2.6) for us, us in terms of U1, U2, Us. 
(ii) We substitute these values in (2.9), and so obtain A, B, C, D in terms of us, U2, Us. 
For future reference, let us solve (2.9) for 73, 74, Ts: 


T; = — DC, T, =C, T, = — AC. (2.11) 
3. An identity satisfied by the Herzberger coefficients for a single surface. Con- 
sider a surface of revolution 


z= f(r), r2 = x? + y?, (3.1) 
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For refraction or reflection at this surface, the angle-characteristic is‘ (if we take the 


origins O, O’ coincident) 


T= (§— (at (a-—ay + — 2s, (3.2) 
from which x, y, z are to be eliminated by the relations 
E-€ dz x n— 7 y 
a cig Cinid tbo tas f(r) 9 scien = f(r) —'? (3.3) 
=a Ox r — r 


It is clear that T will be a function of the two quantities 
o=3[E-&)+-n)*]), w=o-s, 


or, in the notation of (2.2), 
W = 0(n? — 2u3)'/? — O(n’? — 2us)'/?. (3.4) 


od = U3 — U4 + U5, 
Here n, n’ are the refractive indices of the initial and final media, and 6, 6’ are +1; 
for refraction we have 06’ =1, and for reflection 06’ = —1. If we take refraction with 


the rays proceeding in the positive sense, we have 


2 
II 
= 
il 
a 
~~. 
wm 
= 


6 = 1, = —1, n= n'. (3.6) 


By (3.4) we have 
T3 = Ty — TyO(n? — 2us)—”, 
Tr = -Ts, (3.7) 
Ts = To + T,0'(n" — 2us)-!”. 


Hence 
r T (n’? — 2u5)1/? 
eas eng Se ters (3.8) 
T4 + Ts (n? = 2u3) 1/2 
where 
k = 1 for refraction, 
(3.9) 
k = — 1 for reflection. 


Let us substitute from (2.11) in (3.8); this gives 








D-1 (n’? — 2u5)}/? 
a-1 we (3.10) 
A-1 (n? — 2u3)'/? 
and so 
D— 1\? 
Us = 3n.? — 3(n? — 2us) a : 1 . (3.11) 


When we substitute this value in (2.6), and at the same time substitute for 73, 74 


from (2.11), we get 





4 Synge, op. cit., 33. 
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D — 1\? 
n’? — (n? — 2us) (——) = 2(C*u, + CDuz + Dus). (3.12) 
To sum up: For refraction or reflection at a surface of revolution, the coefficients A, C, D 
are connected by the identity (3.12). 
If B=0, then A = D—, and (3.12) simplifies to 


n’? — n®D? = 2C(Cu,; + Duy). (3.13) 

As an alternative procedure we may use the fact that T is of the form 
iad T= (5 — F(x), (3.14) 
x= [(E- 8)? + @— a) I/6 - £')% (3.15) 


This is evident from (3.2) and (3.3); the form of the function F depends on the form 
of the surface. On differentiating (3.14) we obtain three equations analogous to (3.7), 
but containing F and its derivative on the right hand sides. If we eliminate these 
quantities we obtain (3.8) and hence the identity (3.12). 

4. The Herzberger transformation for a sphere. Let us take the origins O, O’ at 
the center of a sphere of radius |r|. The angle characteristic for refraction or reflection 
at the sphere is 

T=+|r|[€-%)?+@—9)? +6 -—%)?)*. (4.1) 
If we suppose the rays incident in the positive sense, all ambiguities of sign are re- 
moved by writing 
T = r(kn’ — n)p', (4.2) 
where 
2 9 9/..19 9 
p = 1 + ———— [kann’ — uy — k(n? — 2us)'/2(n’"? — 2us)/?]. (4.3) 
(kn’ — n)? 
Here r is positive if the rays are incident on the convex side, and negative if they are 
incident on the concave side; k=1 for refraction and k = —1 for reflection. All roots 
are positive. 

We have then 

(n’? — 2u5)!/? 
, 
(n? — 2us)!/? 
T, = — r(kn’ — n)—'p7!/?, (4.4) 
(n? — 2u3)'/? 


(n’? — 2us)'/? 








T; = kr(kn’ — n)—"'p—/? 








Ts = kr(kn’ — n)—'p—"/? 


It is evident that 


i = Ti, (4.5) 


and so, by (2.9) and (2.10), 
B=0, A = 2D", (4.6) 


5 Synge, op. cit., p. 36. 
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We now solve (4.4) for us, us in terms of 73, 74, obtaining 


> Ts 
us = k(n? — 2u3) - 


ps 2 
T. T? 


9 
7? 
—+—-+ h(n’? +n»), 


Ts (4.7) 
us = $n’? — 3(n? — 2us) —- 
TY 
Substitution of these values into (2.5), (2.6) gives 
u2T, = 3(n? + n’) 77 + n?T3T, — 37’, 
I 1 porn? 1 9 2 (4.8) 
Uy = UeT3 + gn'?*T4 — 5n*T73. 
These are two equations for 73, 74; they may be written 
Te= OH 2[ a2 — 3(n? + n”)T, LT |, (4.9) 
Ti(n?— n”)?+ 47; [p?— 37(n?+ n?)| + = 0, (4.10) 
where (in Herzberger’s notation) 
p? = 2n?u, — us. (4.11) 


Solving (4.10) we get, after some simple reductions, 
C = Tp = r[60,(n? — p2/r?)/2 + 04(n’? — p?/r?)1/2], (4.12) 


where 9; and 62 are each +1, for the moment undetermined. We remove the ambiguity 


of sign by considering the case £=7=0, so that by (2.8) &’=Cx, n’ =Cy. It is evident 
from elementary considerations that C has the same sign as (n—kn’)/r. Therefore 


6,=1, 6.= —k, and so in general 


C = r[(n? — p?/r?)1/2 — k(n’? — p2/r?) 2], (4.13) 
By (2.11) and (4.9) we have 
D = n-*[p?/r? + k(n? — p2/r?)/2(n’2 — p?/r?)'/2 — wC]. (4.14) 


We verify that if x=y=0, then D=kn'/n, as it must be by (2.8) from elementary 
considerations. It is easy to check that (3.13) is satisfied by (4.13), (4.14). 
For the case of refraction (k=1) the formula (4.13) agrees with Herzberger’s 


equation (36)!, except for a reversal of sign. 
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ON THE FORCE AND MOMENT ACTING ON A 
BODY IN SHEAR FLOW* 


By YUNG-HUAI KUO (California Institute of Technology) 


Recently, H. S. Tsien solved the problem! of a Joukowsky airfoil in a steady, two- 
dimensional flow of constant vorticity distribution. It is interesting to note that the 
hydrodynamical forces can be expressed in a form similar to the well known Blasius’ 
theorem, involving contour integration of the complex potential function. The follow- 
ing derivation of the formulae is believed to be simpler than that of Tsien. 

1. Equations of motion. Let u and v be the velocity components parallel to the 
x- and y-axis, respectively. In the case of two-dimensional steady motion, the Eulerian 
dynamical equations are: 


Ou dv ov «(Ou 1 dp 
u—+v—-—vi—- — =— (1.1) 
Ox Ox Ox dy p Ox 
Ou Ov do Ou 1 dp . 
e—+9—+s8(——-—}) 2 -——») (1.2) 
oy oy Ox dy p oy 
where ? is the pressure and p, the density of the fluid. The equation of continuity is 
Ou dv 
an fp me oe @ (1.3) 
Ox oy 


For the type of shear flow considered by Tsien,' the vorticity is constant every- 
where in the field and equal to —k. Thus 


Ov Ou 
——-—=-k k>0O. (1.4) 


At the first sight, it seems that the problem might not be definite as one has four 
equations for three variables. By eliminating p between Eqs. (1.1) and (1.2), however, 
the result can be reduced to Eq. (1.3) by means of Eq. (1.4). This shows that any solu- 
tion which satisfies Eqs. (1.3) and (1.4) is consistent with Eqs. (1.1) and (1.2). 

To simplify the problem, the solution is written in the following form: 


u=ky+u', (1.5) 
v=. (1.6) 
Then Eqs. (1.3) and (1.4) reduce to 

ou’ dv’ 
— + — = 0, (1.7) 
Ox oy 

ov =—s us’ 
—--—=0. (1.8) 
Ox oy 


* Received June 21, 1943. 
1H. S. Tsien, Symmetrical Joukowski airfoils in shear flow, Quarterly Appl. Math., 1, 129 (1943). 
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These equations are satisfied by 


F) ay 
wa y= — —; (1.9) 
oy Ox 
or 
deg dy 
Pan A jmu: (1.10) 
Ox oy 
where y¥ and ¢ are the imaginary and real parts of the complex potential F(z); namely, 
¢ + wy = F(z), 2= 2+ ty; (1.11) 
and 


u’ — iv’ = w'(z). (1.12) 


For a given problem the function F(z) is so determined that the velocity component 
norimal to the contour of the body is zero. 
3y virtue of Eqs. (1.4), (1.5), and (1.6), Eqs. (1.1) and (1.2) give 
Pp 49 LU 
p= —— a? — phu'y + pkey, (1.13) 
where g” =u’?+-y”, and the constant of integration is absorbed in y. 


2. Force and moment. If the motion is two-dimensional and steady, the compo- 
nents of the hydrodynamical force and moment? acting on the body are given by 


X= -$ pdy - of u(udy — vdx), 2.1) 
g riztog v(vdx — udy), (2.2) 


M -$ p(xdx + ydy) — p (— v’xdx — u*ydy + uvydx + urxdy), (2.3) 


Y 


where the contour integrals are taken along a closed curve containing the body. Using 
Eqs. (1.5), (1.6) and (1.13), the above equations can be written as: 


p ' 
= og [(u’? — v'*)dy — 2u'v'dx| — oh [(y + u’y)dy —v’ydx], (2.4) 


p . ] 
Y=- 2g [(u!? — o'2)dx + 2u'v'dy| + kh [(y — u’y)dx —v'ydy|, (2.5) 
p 
M:= — Re Ba suds] 
+ oh [(w — u'y)(xdx + ydy) — (v’'yx — 2u’y)dy + v’ y*dx]. (2.6) 


*W. F. Durand, Aerodynamic theory, vol. 2, Springer, Berlin, 1935, pp. 31-33. 
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If only bodies with closed boundary are considered, no sources can exist within the 
field of flow. Then the stream function y is single-valued, and 


¢ ydx -$ x(v’'dx — u’dy), 
¢ ydy -$ y(v'dx — u'dy). 


From these relations, it is not difficult to deduce 


X=- -§ [(u’? — v'2)dy — 2u'v'dx], (2.7) 


] 


I 
+ 


“ = ¢ [(u!? — 9'%)dx + 2u'v'dy] 


+ pk [v'(xdx — ydy) — u'(ydx + xdy)], (2.8) 


p 
— Re - - sw'ds| 
2 


pk 
2 


M 


H 


+2 $ [-w'l(a— yay + 2ayde} + of at — yds — 2eydy}]. 02.9) 


These at once suggest the following alternative expressions: 


i 
X-iW¥= = $ w'dz + i Im E $ wads], (2.10) 
» ike\? 
M = - Re[ =f (w -*) as]. (2.11) 


Eqs. (2.10) and (2.11) may be regarded as an extension of Blasius’ theorem. They 
can be easily identified with the expressions given by Tsien.! The calculation of force 
and moment, however, can be simplified to a certain extent by using these new ex- 


pressions. 
The writer wishes to thank Dr. H. S. Tsien for the use of his paper before publica- 


tion and for his helpful discussions. 


and 
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A CHART FOR PLOTTING RELATIONS BETWEEN 
VARIABLES OVER THEIR ENTIRE REAL RANGE* 


By L. H. DONNELL (Illinois Institute of Technology) 


The following simple method of graphical respresentation, covering the entire real 
range, seems rather obvious and the writer has found it useful for many years; how- 
ever he has never seen it described in the literature. It consists of an ordinary Car- 














0 0 - / 
r 
f 
Y 
/ pf / 
yo HO 
of -/ 
/ 
a 
$ Xx 
om * - ° / P 
x 
Fic. 1. 


tesian plot over the range —1 to +1 of each variable, with adjoining Cartesian plots 
of the reciprocals of the variables from 0 to —1 and from +1 to 0, arranged as shown 
in Fig. 1. This is evidently equivalent to a p!ot over the range — © to + for each 
variable. It is also evident that curves and their slopes will be continuous over the 
dividing lines between the two kinds of plots if the functions represented and their 
first derivatives are continuous at these points. 





* Received March 22, 1943. 
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To illustrate the method, the power relation x =y* has been plotted in Fig. 2 for 
values of x and y between — © and +, and for various values of n. 



































Fic. 2. Plot of x=y". 


In many applications only a part of such a chart would be required. For instance 
only the upper right hand quarter would be needed to cover all positive real values 
of the variables. 

In some cases a change of variable will make the resulting plot much more useful. 
For instance a plot of the relation x = 10*y*, (see Fig. 1) follows the coordinate axes 
closely and would be nearly indistinguishable from many other functions, such as 
x =10'y® or x =10!*y5. On the other hand, a plot of x=y”, where y’ =100y, gives a 
much more illuminating and characteristic picture of the function. Any change of 
variable which makes the resulting curve pass through or near such points as (1, 1), 
(1, —1), (—1, 1) or (—1, —1) would accomplish this purpose in most cases. 
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THE LINES OF PRINCIPAL STRESS IN THE PLANE 
PROBLEM OF PLASTICITY* 


By W. S. AMENT (Brown University) 


Consider a state of plane strain in an incompressible plastic body yielding under a 
constant maximum shearing stress. }. Boussinesq' has shown that the lines of prin- 
cipal stress then form an “equiareal pattern,” i.e. from the two families of lines of 
principal stress individuals can be selected so as to render equal in area the meshes 
formed by these lines. In a recent paper? M. A. Sadowsky has stressed the importance 
of this result and has coined the term “equiareal pattern.” The present note aims at 
establishing the relation between Boussinesq’s result and a theorem concerning the 
lines of principal curvature on certain Weingarten surfaces. 

If the lines of curvature are chosen as parametric curves and « and x’ denote the 
principal curvatures corresponding to the directions of »=const. and u=const. re- 
spectively, the Mainardi-Codazzi relations take the form’ 





0 (log G) Ox’ 
=— (10 = ——) 
ou k—«’ du 
(1) 
rs] (log E) 2 Ox 
a Oo i ee 
Ov k—«’ Ov 


Consider now the Weingarten surfaces for which the difference of the principal curva- 
tures has a constant value. Elimination of x and x’ between the equations (1) then 


leads to 
2 


(log EG) = 0. 





Oudn 


Hence EG =f(u)g(v). A transformation of the type #=4a@(u), 1=0(v) only relabels the 
parametric curves but does not affect their geometric properties. Define @ and i by 


da oa do _ = 
_* Vf(u), — g(v). 


For these new parameters EG=1, i.e. the meshes formed by the parametric curves 
corresponding to two sets of equidistant values of # and @ are equal in area. The lines 
of principal curvature on.a Weingarten surface with x—x'=const. therefore form an 
equtareal pattern. 

The relation between this theorem and Boussinesq’s result is immediate. Intro- 
duce rectangular Cartesian coordinates O, x, y, 2, the plane O, x, y having the orienta- 
tion of the plane of strain. The normal stresses o,, a, and the shearing stress 7 then 
can be derived from a stress function F according to 


* Received Nov. 6, 1942. 

1C, R. Ac. Sci. Paris, 74, 242 (1872). 

2 Trans. Am. Soc. Mech. Eng. 63, A-74 (1941). 

3 See for instance C. E. Weatherburn, Differential geometry, vol. I1, Cambridge 1930, p. 52. 





W. S. AMENT 
oF °F oF 
z =_l oO , Oy = _ , - = _ . 
ay? ax? axdy 


The yield condition 
(o, — oy)* + 41? = const. 


furnishes the following differential equation for F: 


OF O2F\2 O°F \2 
(— - =) +4/( ) = const. 
0x" oy" Oxdy 
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(2) 


Consider the surface S defined by 2(x, y)=aF(x, y), where a is an arbitrary small 
constant rendering small the slope which the tangential planes of S have with respect 
to the plane O, x, y. If terms of the second order in 0F/dx and 0F/0y are neglected, 
the difference of the principal curvatures of S is seen to equal the square root of the 
left hand side of (2) multiplied by a. The surface S therefore is a Weingarten surface 
of the type considered above. The orthogonal projections of its lines of principal 


curvature on the plane O, x, y are the lines of principal stress. 





CONSERVATION OF SCHOLARLY JOURNALS 


The American Library Association created in 1941 the Committee on Aid to Li- 
braries in War Areas, headed by John R. Russell, the Librarian of the University of 
Rochester. The Committee is faced with numerous serious problems and hopes that 
American scholars and scientists will be of considerable aid in the solution of one of 


these problems. 


One of the most difficult tasks in library reconstruction after the first World War 
was that of completing foreign institutional sets of American scholarly, scientific, and 
technical periodicals. The attempt to avoid a duplication of that situation is now the 
concern of the Committee. 


Many sets of journals will be broken by the financial inability of the institutions 
to renew subscriptions. As far as possible they will be completed from a stock of 
periodicals being purchased by the Committee. Many more will have been broken 
through mail difficulties and loss of shipments, while still other sets will have disap- 
peared in the destruction of libraries. The size of the eventual demand is impossible 
to estimate, but requests, received by the Committee already give evidence that it 


will be enormous. 


With an imminent paper shortage attempts are being made to collect old periodi- 


cals for pulp. Fearing this possible reduction in the already limited supply of scholarly 
and scientific journals, the Committee hopes to enlist the cooperation of subscribers 
to this journal in preventing the sacrifice of this type of material to the pulp demand. 
It is scarcely necessary to mention the appreciation of foreign institutions and scholars 


for this activity. 


Questions concerning the project or concerning the Committee’s interest in par- 
ticular periodicals should be directed to Dorothy J. Comins, Executive Assistant to 
the Committee on Aid to Libraries in War Areas, Library of Congress Annex, Study 
251, Washington, 25, D. C. 

AMERICAN LIBRARY ASSOCIATION, 
Committee on Aid to Libraries in War Areas. 














